Correspondence principle for many-body scars in ultracold Rydberg atoms

The theory of quantum scarring -- a remarkable violation of quantum unique ergodicity -- rests on two complementary pillars: the existence of unstable classical periodic orbits and the so-called quasimodes, i.e., the non-ergodic states that strongly overlap with a small number of the system's eigenstates. Recently, interest in quantum scars has been revived in a many-body setting of Rydberg atom chains. While previous theoretical works have identified periodic orbits for such systems using time-dependent variational principle (TDVP), the link between periodic orbits and quasimodes has been missing. Here we provide a conceptually simple analytic construction of quasimodes for the non-integrable Rydberg atom model, and prove that they arise from a"requantisation"of previously established periodic orbits when quantum fluctuations are restored to all orders. Our results shed light on the TDVP classical system simultaneously playing the role of both the mean-field approximation and the system's classical limit, thus allowing to establish a rigorous link between the wavefunction scarring in the Rydberg atom chains and the single-particle quantum systems.


I. Introduction
Quantum scars provide a surprising connection between single-particle quantum billiards and their classical counterpart. 1,2 The understanding of quantum scars is built upon the existence of unstable classical periodic orbits and atypical eigenstates. The former was first observed by studying the spreading of wave packets initialised along such orbits 1,3 (see also Ref. 4). In the limit → 0, the classical orbits are recovered, as stipulated by the correspondence principle. 5 Furthermore, some atypical "scarred" eigenstates were found to concentrate around such orbits, in contrast to typical eigenstates at a similar energy, which tend to be well-spread across the whole configuration space. In the fully quantum case, the proof for the existence of scarred eigenstates relied on approximate eigenstates or "quasimodes", which are constructed to be strongly localised around a classical periodic orbit. For example, in the celebrated Bunimovich stadium problem, 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. 2,6 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 atypical properties. [6][7][8] Recently, similar non-thermal states have been identified in interacting quantum systems. These states, named quantum many-body scars by analogy, have been linked to long-lived oscillations measured in an experiment on a chain of Rydberg atoms. 9-16 Since then, an increasing variety of models have been associated with various aspects of many-body scarring.  In this work we focus on the "PXP" model -an idealised effective model of the Rydberg atom experiment 49,50 -which can be formally expressed in terms of interacting spin-1/2 degrees of freedom, without an "obvious" semiclassical limit. Inspired by the theory of single-particle scars, an immediate question arises: what are the periodic orbits and quasimodes in the PXP model?
Unlike single-particle billiards, the equivalent of the classical trajectory in the many-body case is less transparent. One proposal makes use of the Time-Dependent Variational Principle (TDVP), 51 wherein many-body dynamics is projected onto a restricted manifold of matrix product states. 52-54 Its application to the PXP model successfully captures the revivals following quenches from specific product states. 11,24 It is thus 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 correspondence principle.
On the other hand, a family of quasimodes for the PXP model have also been independently constructed by the so-called Forward Scattering Approximation (FSA). 10,55 This method makes it possible to extract properties of scarred eigenstates in systems much larger than those accessible by exact diagonalisation (ED). However, this way of building quasimodes presents some deficiencies, both on technical and conceptual levels. Indeed, while the computational cost of generating the FSA states is better than the exponential growth of the Hilbert space, it is still polynomial in the system size. Furthermore, the link between the FSA and the TDVP is obscure, thus making it difficult to relate the quasimodes with the classical limit of the model, in a manner that had been achieved for single-particle scars.
In this paper we introduce a new way to form quasimodes which resolves the aforementioned problems. Our construction relies on building a subspace symmetrised over permutations within each of the sublattices comprising even and odd sites in a chain. This approach can be thought of as a "mean-field" approximation, and makes it possible to obtain closed-form expressions for the projections of states and operators in this symmetric subspace. This enables us to numerically approximate highly-excited eigenstates of chains of length N 800, far beyond other methods. Further, we show that the path-integral quantisation of the TDVP topological space with a proper resolution of the identity exactly reproduces the symmetric subspace in the thermodynamic limit. While adding quantum fluctuations to the TDVP has been approximated to first order, 56 our method includes quantum fluctuations to all orders. This correspondence highlights the intriguing fact that TDVP in the PXP model serves as both a mean-field approximation and a classical limit taken together. We also show that the scarred quasimodes are localised around the TDVP periodic orbit, in a way that is reminiscent of single-particle quantum scars. Finally, performing time evolution in the symmetric subspace strongly suggests that the classical periodic orbit found in TDVP is stable to the addition of quantum fluctuations in the thermodynamic limit, at least to the "mean-field" level as described above.

II. Rydberg blockade and quantum scars
The PXP model is an effective model of a chain of atoms in the Rydberg blockade. 49 Denoting by |• and |• the ground and excited states of each atom, the Hamiltonian is where X n = |• •|+|• •| is the Pauli X and P n = |• •| is a local projector. Here N denotes the number of sites in the chain and, unless specified otherwise, we assume periodic boundary conditions (PBCs), N + 1 ≡ 1. The projectors make the model in Eq. (1) constrained: a site can be excited only if both adjacent sites are in the ground state, e.g., the transition |. . . ••• . . . → |. . . ••• . . . is forbidden. This causes a fracturing of the Hilbert space, which has a strong impact on its dynamical properties. 57-62 Below we focus on the largest connected sector of the Hilbert space which contains no state with adjacent excitations. The PXP model in Eq. (1) is non-integrable and thermalising, 10 but it exhibits periodic revivals when quenched from certain initial states such as the Néel state, |Z 2 = |•••• . . . . The observed revivals in the density of domain walls 9 are remarkable given that the |Z 2 state effectively forms an "infinite-temperature" ensemble for this system, for which the "strong" Eigenstate Thermalisation Hypothesis (ETH) predicts fast equilibration. 63-66 Unlike many-body localised 67-69 and integrable systems 70 which strongly break ergodicity, the PXP model breaks it only weakly. 71,72 The breaking of ergodicity in the PXP model occurs most prominently due to the existence of a band of O(N ) scarred eigenstates, which have equal energy separation, anomalously high overlap with |Z 2 state, and lowest entanglement among all eigenstates. 55 These states are the representatives of scarred towers of eigenstates (as can be seen in Fig. 2(a) below), which contain additional atypical eigenstates that exhibit clustering around the same energies and whose finite-size scaling of inverse participation ratio is different from that of the thermalising eigenstates. 10 The scarring properties in all of these special eigenstates are explained by the existence of quasimodes: simpler states with an underlying structure -the analogues of "standing waves" in billiard systems -which represent concentration points of the special eigenstates.

III. Construction of quasimodes
We construct quasimodes for the PXP model by first defining a symmetric subspace K from 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 which "shuffle" the sites in each sublattice. This approach is somewhat reminiscent of the fully symmetric subspaces in the context of quantum de Finetti theorem 73 and in studies of fully-connected models. 74,75 However, the key difference here is that the permutation symmetry is broken to the sublattice level, and many of the permutation shuffles violate the Rydberg constraint and therefore have to be excluded, which makes the analysis non-trivial. An illustration of the construction of K for a PXP model of size N =8 is presented in Fig. 1.
We define an orthonormal basis for K from the symmetric combinations of members of each class (n 1 , n 2 ), |(n 1 , n 2 ) = 1 #(n 1 , n 2 ) where #(n 1 , n 2 ) is the size of the equivalence class.
The class sizes can be calculated analytically because the configurations can be generated by a recursive procedure or transfer-matrix problem which 'glues' sites (in the basis) onto the boundary of a smaller version of the same counting problem. Because the constraint on gluing an additional • site is only sensitive to the state of the site immediately on the boundary, we can erase the information about the configurations further away and represent them by a count of all those configurations which are equivalent on the boundary. This produces a method for calculating #(n 1 , n 2 ) using dynamic programming in only polynomial time. It can be shown 76 that the class sizes admit the closed-form expression in the case of PBC, with an analogous expression for open boundary condition. For simplicity, we assumed N is even. Using Hermiticity and the sublattice exchange symmetry, all matrix elements of H in the symmetric subspace can be found by examining only the (H) (n1−1,n2) (n1,n2) = (n 1 −1, n 2 )| H |(n 1 , n 2 ) matrix elements. These can be found using the idea that for each configuration in (n 1 , n 2 ) there are precisely n 1 ways to remove an excitation from sublattice 1. Accounting for normalisation and using the closed form Eq. (3) yields a remarkably simple effective Hamiltonian within K:

IV. Properties of quasimodes
The quasimodes, or approximations to eigenstates, are formed by diagonalising the PXP Hamiltonian projected to K, Eq. (4). Despite the apparent simplicity of the matrix elements, we are unaware of an analytical method to diagonalise Eq. (4) and we perform this step numerically. In Fig. 2(a) we compare the quasimodes against the exact eigenstates for system size N =32. For both the eigenstates and quasimodes we plot their mean energy and overlap with |Z 2 state. N + 1 scarred quasimodes, indicated with red crosses, stand out by closely matching the scarred eigenstates, which are responsible for the observed quantum revivals. These scarred quasimodes have previously been constructed via the FSA. 10 In contrast, the present method yields a larger set of O(N 2 ) quasimodes. The remaining quasimodes, indicated with the blue pluses in Fig. 2(a), also follow a regular pattern.
We measure the localisation of the quasimodes among energy eigenstates using participation ratio, PR = E E|φ a 4 , which represents the inverse of the number of eigenstates with which a state has 'typical' overlap. For random states this is 1/D, with D the Hilbert space dimension. In Fig. 2 of all the quasimodes. The N +1 scarred quasimodes have very large PR compared to a typical random state, with appreciable overlap with only a small number of eigenstates. The remaining quasimodes are also significantly more localised than typical states. In Fig. 2(c) we show the scaling of PR with N , focusing on the scarred quasimode closest to the middle of the spectrum (E = 0) in the zero-momentum sector, alongside nonscarred quasimodes with energies 0<E< √ N . This avoids the edges of the spectrum where the density of states is much smaller. From this we conclude that the scarred quasimodes are localised over a small number of eigenstates compared to the density of states.
Furthermore, the simple structure of the symmetric subspace allows us to characterise the errors in our eigenstate approximations beyond the limits of ED. We can calculate the energy variance matrix in the symmetric subspace, σ 2 H =K(H 2 )−(K(H)) 2 , whereK is a superoperator which projects linear operators on the full Hilbert space down to the symmetric subspace K. It can be shown that the energy variance matrix is diagonal in the −20 −10 0 log 10 | Z 2 |E i | 2 Scarred quasimodes (red crosses) have low energy variance density compared to typical energy scales, indicating that the dynamics restricted to K remains close to the exact evolution to relatively long times. (b) Energy variance in the symmetric basis, Eq. (2), as a function of scaled n1 and n2. Superimposed (in red) is the probability distribution on n1 and n2 for a scarred quasimode in the middle of the spectrum (with E ≈ 27.28).
symmetric basis and has the following closed form: 76 wheren 1 ,n 2 are the operators that count excitations on their respective sublattices, e.g.,n i |(n 1 , n 2 ) = n i |(n 1 , n 2 ) for i = 1, 2. As expected, among the quasimodes, we find that the scarred ones are those with the smallest energy variance, especially those that have high overlap with |Z 2 state. This can be seen in Fig. 3(a), which slows the energy variance for individual quasimodes in system size N =140. The variance is linear in N , representing constant density of energy fluctuations, but for the scarred quasimodes it is much smaller than the typical energy scales. The average energy fluctuations in a quasimode is asymptotically given by 1/9, while the same quantity restricted to scarred quasimodes is found numerically to be ≈ 0.0049. This suggests that the time evolution approximated within this subspace will be reasonably accurate for local quantities such as local observable expectations for relatively long times. The typical scarred quasimode support can be visualised on an (n 1 ,n 2 ) triangle -see Fig. 3(b), which reveals that in these cases the quasimodes tend to "avoid" the region with n 1 and n 2 both large where the energy variance is high. This is analogous to the classical trajectory in the TDVP avoiding "high-leakage" regions of the phase space. 11

V. Quantum dynamics
Our construction of the quasimodes can also be used to approximate the dynamics of the system. This is illus-trated in Fig. 4, which displays a comparison with ED as well as variational (TDVP) results, for a quench from the Z 2 Rydberg crystal. TDVP approach is based on a restricted class M of matrix product states |Ψ(θ, φ) with bond dimension equal to 2, which can be parametrised using matrices 9,11 where i labels the lattice sites or, in an infinite system, sites within a unit cell (for the |Z 2 quench, i = 1, 2). The states are asymptotically normalised, Ψ (θ, φ)|Ψ(θ, φ) → 1, as N → ∞. 11 The time evolution of the TDVP angles in Fig. 4 was obtained using the equations of motions for the infinite system. 24,76 The differences between the dynamics within K and the TDVP dynamics on M can be seen clearly in Fig. 4(a), which shows the von Neumann entropy for an equal bipartition of the system. This quantity, although nonlocal, can also be efficiently calculated using the present method. 76 We observe that entropy captures the differences between approximation schemes well before the relaxation time set by the energy variance. In contrast, local quantities such as local expectation values are very consistent -see Fig. 4(b) which shows the density of Rydberg excitations in the system, n/N . The accuracy of these different methods can be assessed using co-moving fidelity densities between pairs of evolutions, such as 1 N log | ψ ED (t)|ψ K (t) | 2 . Both the TDVP and symmetric subspace approaches generate a very low error density. The dynamics in the symmetric subspace has been simulated in large system sizes (for N ≤ 720 76 ), and the results are consistent with the revivals of the Néel states becoming perfect in the thermodynamic limit. However, there is a slight deviation from the TDVP orbit, 11 in that there is no perfect state transfer between |•••• . . . and the Néel state with the two sublattices flipped, |•••• . . . .
The co-moving fidelities in Fig. 4(c) underscore there is a curious similarity between the evolutions within TDVP and symmetric subspace K. The variational class of states considered in the TDVP approach can be viewed as a Gutzwiller projection of spin coherent states with a unit cell of two sites, where the projector eliminates nearest neighbour excitations. 9,11,24 The states in this subspace M are parametrised by two angles which represent probabilities for a site of each sublattice to hold an excitation. Recall that the basis for K consisted 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 N . 76 This implies that the dynamics within K has a leakage rate, or instantaneous error, which is bounded above by that of the TDVP evolution. Next we will investigate this connection in more detail and show how K and M are related by quantisation.  The entropy growth shows that for K, the initial state does not recur exactly leaving some residual entropy compared to the more weakly entangled dynamics on M. (b) Local observables following the quench are well reproduced. Both TDVP and K subspace approximations fail to capture a slow decay in the oscillations for t 5. (c) Pairwise comparisons between the exact dynamics, projected dynamics on K and TDVP dynamics on M. This illustrates that the dynamics on K is always more accurate than TDVP on M.

VI. Quantisation of periodic orbits and the correspondence principle
Our final step is to establish a correspondence between K and TDVP. The latter is typically conceived as a variational approximation method -how do we reinterpret it as a classical limit? Conversely, how can we best quantise this classical system, given that quantisation is often a non-unique recipe? These questions are conveniently addressed in a path integral framework. Thus, we first carefully formulate the path integral over our constrained space before arriving at the correspondence principle that relates M and K.

A. The variational principle and dequantisation
A path integral over a classical system is usually used to define a quantum propagator which associates amplitudes to pairs of path endpoints. This propagator may be viewed as a linear operator over the span of the classical space. However, in the present case, we follow a reverse procedure: we use the path-integral formulation to "dequantise" a system by writing its propagator as a path integral and then take its constituent classical state space and action functional as defining a classical system. When the TDVP state space furnishes a resolution of the identity, the TDVP retraction is equivalent to such a dequantisation. 76 This process provides a correspondence between a classical system such as M and a quantum system in its linear span -for M, this is the Rayleigh-Ritz projection of the dynamics in the symmetric subspace K. This is clearly the most natural outcome, as it meets our expectation that the different representations of quantum mechanics -i.e. position, momentum or coherent-state -are precisely that: representations of the same object. If the same variational principle is taken for restricting to a vector subspace, such as our symmetric subspace K, then the Rayleigh-Ritz method is exactly recovered. It is natural that both the quantised and dequantised system satisfy this same same kind of variational optimality.
Crucially, the prescription above relies on being able to associate a measure µ to some family of paths. This can be achieved by equipping the classical state space with an integration measure, forming a structure we will refer to as a frame. 77 This structure is characterised by the frame operator, using the parametrisation of states in Eq. (6). For (M, µ) to constitute a genuine frame, S µ must be bounded and positive definite. It is this operator which must reproduce the identity to have a well-defined quantisation. Unfortunately, our state space M does not come with a natural choice of µ, unlike many of the familiar manifolds of quantum states, such as spin coherent states and matrix product states, 78 where µ is simply given by the Haar measure for a transitive group action. 79 Even a measure which works for the space of constrained coherent states where every site has independent variables (θ j , φ j ), found in Ref. 11, does not survive when these are restricted to two-site periodicity.
The strategy we will follow is to start from a measure that produces a frame over M which is at least reasonably well-behaved, although it may not resolve the identity. The corresponding frame operator S µ will turn out to be bounded and positive definite for all finite N , albeit merely positive semidefinite in the limit N → ∞. This can then be used to construct a transformed state space, M, that does resolve the identity, and thus quantise naturally. Formally, we construct this new frame (M,µ) from our existing one by replacing the states attached to each point of the space according to the linear action, We refer to this as a frame transformation in analogy to a basis transformation. This process is well defined whenever S µ is bounded and positive definite, and the corresponding frame operator is the identity operator. In this case, the thermodynamic limit can be regularised such that it quantises naturally to a space with only some "non-physical" states removed. The removal of non-physical states is justified because, by construction, M is meant to be effective at describing the Néel state quench, while in other parts of the phase space the leakage rate or geometric error may be large. 11,24 For finite sizes, we will show that the transformed space is a moral equivalent to M, in the sense that they are highly similar around the "good" trajectories, while potentially disagreeing about the "bad" trajectories. Furthermore, we will show the differences between the two shrink as N increases, until in the thermodynamic limit they are equivalent almost everywhere. Our starting point is the gauged coherent state parametrisation for M, Eq. (6), and we want to find a measure such that S µ is bounded. The choice of measure is to an extent arbitrary, since its details will be transformed away later, but the challenge of ensuring that S µ is bounded can be understood from the following heuristic argument. Consider the expectation value of S µ for a point x of M. A contribution to this integral will be significant if it is not too distinguishable from x, and therefore the integral can be estimated by the volume of the "fuzzy" set indistinguishable from x. For almost all x and as N increases, the fuzzy neighbourhood retracts around x to a point, with a typical fuzzy size . This phenomena also occurs on the manifold of spin-coherent states with increasing spin. 80 For x along the edges (i.e. θ 1 = π), however, it retracts to a line, with a typical fuzzy width O(1/ √ N ). The line-like contributions to S µ appear at a higher order in N , due to their different volume scaling, thus spoiling a bounded positive-definite limit. This issue can be avoided by adopting for our measure Relative to the product of spherical measures, the additional cos 2 θi 2 is chosen to vanish rapidly enough at the edges of the phase space that the singularity is removed.
We proceed by examining matrix elements of S µ in the symmetric basis of K. Only the diagonal matrix elements S µ n1,n2 = (n 1 , n 2 )| S µ |(n 1 , n 2 ) are non-zero. For periodic boundary conditions these are where (x) (n) is the rising factorial or Pochhammer polynomial. 81 As mentioned previously, we see that the frame operator S µ is bounded and positive definite for every finite N , but it is merely positive semidefinite in the N →∞ limit. The inverse is only ill conditioned near the diagonal in the (n 1 ,n 2 ) plane, but this subspace dynamically decouples from the rest of K since the Hamiltonian matrix elements involving this subspace vanish asymptotically. Moreover, when represented on M, it concentrates for large N around the corners, (θ 1 , θ 2 ) = (±π, ±π), which are not physically relevant for the Néel state quench, besides being measure zero. We can form a convergent series for the inverse by regularising it with a scalar, (S µ + ) −1/2 , which has a projector onto K with the decoupled space removed for its limit. This allows us to 16 32 64 128 256 N 5 10 20 Figure 5. (a) Fubini-Study distance DFS between original M and transformed M state spaces, averaged over the phasespace with measure µ. (b) Illustration of the integral curves of (red) the frame-transformed dynamics and (grey) the original dynamical system, for N = 128. The flows are highly similar except towards the corners. The black line shows the evolution following a Néel state quench in the transformed system. define the frame-transformed space M even in the thermodynamic limit.
We find that almost all of the phase space is left essentially unaffected by the frame transformation in the thermodynamic limit. 76 This occurs because the frame operator S µ is a function of the densities n 1 /N , n 2 /N , with a positive radius of convergence, except in the thermodynamic limit when approaching the diagonal in the (n 1 , n 2 ) plane. The fluctuations in the density observables decay as O(1/ √ N ) over almost all of M because the correlation length of connected correlation functions only diverges towards the corners where the states are dominated by configurations composed of a small number of large Néel domains. This means that S µ asymptotically acts essentially as a scalar on each of the points of M, 76 and can therefore be absorbed into the measure. Asymptotically, the two state spaces, M and M, are equal almost everywhere, and their convergence is numerically demonstrated in Fig. 5(a), which shows the Fubini-Study metric or Bures distance 82 D FS -which measures the angle between two states as projective rays -averaged over the phase space. In the large-N limit, we find D FS ∼ 1/ √ N , which supports the asymptotic equivalence between M and M. The change in measure is unable to affect the behaviour of the classical dynamical system, which can be formulated entirely independent of the measure.
Wherever the fluctuations are small, the frame operators acts approximately as a scalar with small corrections due to density fluctuations. Hence, even for finite system sizes, the transformed space is very similar to the original space except around regions of the phase-space which are considered non-physical. This is shown in Fig. 5 where the grey lines are the original dynamical system and the red lines are the frame-transformed dynamics. The vector fields are for finite-sized systems and are calculated from the action principal represented in the symmetric basis for K. This method avoids numer-ical instability found in the canonical treatment 51 around singularities in the tangent-basis Gram matrix. We observe in Fig. 5 (b) that the Hamiltonian vector fields agree around the good trajectories and diverge from one another on approaching the corners. The black line is an integral curve for the frame-transformed system starting from the Néel state, thus showing that the trajectory which underlies the quantum scar dynamics is retained.
To summarise, we have shown there is a measure which quantises the original space M to K, in the thermodynamic limit and after removing the physically-irrelevant subspace. This was achieved by constructing a sequence of transformed spaces for finite sizes which has a common N → ∞ limit as the original space. Each finite N transformed space quantises naturally to the Hilbert space K examined in the previous sections, in which many quantities of interest may be calculated efficiently. Even for finite system sizes, the transformed space is very similar to the original space except around regions of the phasespace which are considered non-physical.

B. Physical interpretation of quasimodes
Having provided a quantisation-dequantisation correspondence between the classical TDVP system and the symmetric subspace K, we now turn to the consequences for our interpretation of individual quasimodes, and by extension the exact eigenstates. We will show that our scarred quasimodes can be viewed as standing waves localised around the Néel state trajectory, making contact with the the phenomena which inspired the term quantum scar in quantum billiard systems. In these systems, it was found that some eigenfunctions anomalously localise around isolated periodic orbits, 1 leading to the picture that the classical trajectory is "imprinted" upon the quantum states, giving rise to wavefunction scarring.
First, we will show how the quasimodes can be viewed as wavefunctions. The frames constructed in the previous section, which resolve the identity, allow us to analyse a quantum state |ψ ∈ K through the inner-product map: This mapping allows us to form a wavefunction representation, ψ(θ, φ), and, conversely, to recover the quantum state from that representation. This serves as a phasespace analogue to a wavefunction, which in the usual Schrödinger picture would be defined only on a Lagrangian submanifold such as configuration space. Squaring this wavefunction creates an analogue to the Husimi Q-representation which is a quasi-probability distribution rather than a true probability distribution because its events are not disjoint. 80 We restrict attention to a Lagrangian section with φ a = 0, as this typically loses no information and the twodimensional space is better suited for visualisation. 24 The Figure 6. Viewing a selected (a) scarred and (b) non-scarred quasimode as wavefunctions over θ1 and θ2 for N =128. Colour scale represents the real part of the wavefunction. Scarred quasimodes concentrate around the classical periodic orbits demonstrating quantum scarring. The phase of the wavefuntion winds along the trajectory, and this winding number is related to the approximately equally spaced energies of the quasimodes. The other quasimodes avoid this periodic trajectory, due to orthogonality, and concentrate around the trajectories with the corners of the square for endpoints. real part of the wavefunction representation in Eq. (11) is shown in Fig. 6, for an illustrative example of a scarred quasimode (a) and one of the remaining quasimodes (b), for N =128. The scarred quasimodes are found to concentrate around the Néel classical trajectory in a manner strikingly reminiscent of wavefunction scarring in the quantum billiards problem. The classical trajectory shown by the black line in Fig. 6 is for the thermodynamic limit. In Fig. 6(a) we observe that the wavefunction spreads out towards the boundary and becomes completely delocalised at the edges, θ i = ±π, where the points are indistinguishable in the transverse direction. Following the classical trajectory, the phase of the wavefunction winds, and integrality of this winding number leads to their approximately equal spacing in energy. We recognise this behaviour from quantisation of regular trajectories: 83 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 scarred exact eigenstates given their strong overlap with the quasimodes.
In contrast, the non-scarred quasimodes are typically found around the orbits connecting corners of the manifold -see Fig. 6(b). The corners represent the zero vector in the quantum picture, around which the geometric error is large. This may account for the increased energy variance and energy delocalisation found previously in such quasimodes.
We note a dissimilarity between our result in Fig. 6 and wavefunction scarring in the fact that the phase-space section is entirely regular. It cannot exhibit chaos because a symmetry argument leaves the two-dimensional φ a =0 section dynamically invariant. 11 The non-scarred quasimodes are distinguished from those we have called scarred not by a chaotic nature but because they do not appear to survive the injection into the full Hilbert space. Because of the greater variational errors en-countered along the other trajectories, we suggest that the variational space could be increased with additional parameters in such a way that the bad trajectories loose stability -leading to the emergence of classical chaos. The Néel trajectory ought to remain regular, producing a phase-space which is mixed or perhaps with isolated regular trajectories. 24

VII. Conclusions
We have introduced a symmetric subspace K which allows for a simple construction of quasimodes in the Rydberg atom model. We showed that these quasimodes are well-localised in the energy eigenbasis and have greatly suppressed energy fluctuations, thereby making them excellent approximations to both the exact eigenstates and the quench dynamics. The simplicity of this approach allowed us to obtain closed form results for many important quantities and our method can be easily extended to account for various perturbations to the PXP model that have been considered in the literature. 12,14,55 Moreover, we expect the analogue of subspace K to be particularly useful in studying two-dimensional scarred models, 39,84,85 and in describing dynamics in other types of constrained models, e.g., quantum dimer models, 86 whose ground states have a similar "symmetrised" structure at Rokshar-Kivelson points. 87 We have argued that our subspace K can be interpreted as a "requantisation" of the variational classical limit, 11 when quantum fluctuations are restored to all orders. This reveals the TDVP classical system is formed of two approximations, one of which is a mean-field type approximation and another which is the classical limit. Using finite size scaling, we showed that the wavefunction revivals survive in the thermodynamic limit in the mean-field regime, but the state transfer between the two Néel state is no longer perfect. This means that the TDVP periodic trajectories are robust to the addition of quantum fluctuations.
These results allows us to complete the analogy between the wavefunction scarring in the PXP model and the single-particle quantum systems. The one outstanding challenge for making contact between the quasimodes and rigorous results on exact eigenstates is the exponential many-body density of states. The subspace constructed here is naturally related to a large approximate symmetry action which shuffles the lattice sites, within the even and odd sublattices -this could provide a way to control the density of states in order to rigorously prove the existence of quantum non-unique ergodicity in the thermodynamic limit of the Rydberg atom chain.

VIII. Acknowledgements
We thank Alexios Michailidis for useful discussions. We acknowledge support by EPSRC grants EP/R020612/1, EP/R513258/1 and EP/M50807X/1. Statement of compliance with EPSRC policy framework on research data: This publication is theoretical work that does not require supporting research data. This work benefited from participation at KITP Follow-Up programme, supported by the National Science Foundation under Grant No. NSF PHY-1748958. In this supplementary material, we present analytical derivations of class sizes and operator matrix elements projected to the symmetric subspace K. We elucidate in detail the connection between K and the TDVP subspace M, proving that K = spanM. We discuss the quality of revivals and state transfer in large systems within the symmetric subspace approximation. Finally, we provide details on the path-integral quantisation in the constrained space.

SI. Analytical derivation of class sizes and matrix elements A. Class sizes
In the symmetric subspace, each basis state is a symmetric superposition of product states in the full (constrained) Hilbert space. All these states have the same number of excitations on each sublattice, and so they belong to the equivalence class (n 1 , n 2 ) N . The basis states of the subspace are then defined as where #(n 1 , n 2 ) N is the size of the corresponding class. In order to compute the matrix elements of H, we now derive an analytical expression for the class size using a recurrence relation for periodic boundary conditions.
We start with the simple case where n 1 = 0, meaning that there are no excitations on the first sublattice. Then the Rydberg constraints are irrelevant and the class size corresponds to the number of ways to place n 2 identical excitations in N/2 sites. This leads to which is simply a binomial coefficient.
For n 1 = 0, we want to obtain a relation between #(n 1 , n 2 ) N and #(n 1 − 1, n 2 ) N −2 . This can be done by counting the number of product states in the class starting with the pattern |••.... in two different ways. This quantity is denoted by D •• (N, n 1 , n 2 ), and can be computed using the operator S i •• , which is a projector on all states having •• on sites i and i + 1. So for any product state |ψ in the full Hilbert space, This allows us to compute Here we used the fact that the sum over the class elements is invariant under translation of 2 sites. The two sums can then be swapped and the sum over j simply counts the number of excitations on the first sublattice, and thus will give n 1 for every state in this class. The second way of counting D •• (N, n 1 , n 2 ) relies on separating the first two sites from the rest of the chain. Any state |ψ ∈ (n 1 , n 2 ) N with the pattern |•• . . . can be rewritten as |•• ⊗ |φ , where |φ ∈ (n 1 − 1, n 2 ) N −2 with its last site not excited. The fraction of states in (n 1 − 1, n 2 ) N −2 satisfying this condition is given by 1 − n2 N/2−1 = N/2−1−n2 N/2−1 . Accounting for the size of this class leads to Putting together equations Eq. (S4) and Eq. (S5) gives the recurrence relation Using Eq. (S2) as the initial condition, it is straightforward to use this recurrence to obtain the general class size

B. Matrix elements
As the PXP Hamiltonian removes or adds a single excitation, |n 1 , n 2 only has non-zero matrix elements with |n 1 ± 1, n 2 and |n 1 , n 2 ± 1 . Using the Hermiticity of H and the exchange symmetry n 1 n 2 , only one of these S2 elements has to be computed. For simplicity, we choose For each element of (n 1 , n 2 ) N , there are n 1 excitations that can be removed to give a product state in (n 1 − 1, n 2 ) N . Because of the orthonormality of the product states, each of these moves simply gives a contribution of 1. The double sum ends up giving n 1 · #(n 1 , n 2 ) N and the matrix element is This can be rewritten using Eq. (S7) as The case where n 1 = N 2 needs to be treated separately. Because of the kinetic constraint, a sublattice can only be fully occupied if the other one is empty, hence n 2 = 0. Eq. (S9) can then be expanded using Eq. (S2) instead of Eq. (S7) to give which is well behaved for all values of n 1 .

SII. Subspace energy variance
The subspace energy variance is defined as whereK is a super-operator which projects linear operators, defined on the full Hilbert space, down to the symmetric subspace K. The matrix elements of (K(H)) 2 can easily be computed from the results of Sec. SI. To obtain the elements ofK(H 2 ), it is helpful to decompose H as where H ± i adds or removes an excitation on the i-th sublattice. This turns H 2 into a product of 16 terms, and we will show that only 4 of them are relevant. To lighten the notation we define the result for each term as . (S15) As these elements are all computed for a system of fixed size N , the subscript referring to it is dropped. An important simplification comes from the fact that H + 1 |(n 1 , n 2 ) = (n 1 + 1) #(n 1 + 1, n 2 ) #(n 1 , n 2 ) |(n 1 + 1, n 2 ) , (S16) which implies that H + i does not leak out of the symmetric subspace. This is not the case for H − 1 , and if it were it would mean that the subspace energy variance is always zero. The difference between these two operations comes from the Rydberg constraint. Indeed, consider a state |ψ ∈ (n 1 + 1, n 2 ) and remove an excitation. As stated in Sec. SI, this gives n 1 + 1 distinct elements in (n 1 , n 2 ), one per excitation that can be removed. Let us denote the set of these states by G ψ ⊂ (n 1 , n 2 ). Formally this means that G ψ := {|φ ∈ (n 1 , n 2 )| φ| H − 1 |ψ = 0}. When applying H + 1 to |(n 1 , n 2 ) , the n 1 + 1 states in G ψ will contribute to |ψ , and the ones in (n 1 , n 2 ) \ G ψ will not. So H + 1 |(n 1 , n 2 ) gives n 1 + 1 contributions to |ψ . As this is valid for any state |ψ in (n 1 + 1, n 2 ), it means that applying H + 1 to |(n 1 , n 2 ) will give a symmetric superposition of states in (n 1 + 1, n 2 ). Taking into account the normalisation factors gives equation Eq. (S16).
For H − 1 the same procedure is not possible. Indeed, the number of excitations that can be added to a state |ψ ∈ (n 1 , n 2 ) is not uniquely defined. It is actually very easy to find an example of two states in the same class to which we can add a different number of excitations.
The main consequence of Eq. (S16) is that the right actions of H + i andK(H + i ) are the same. This also holds for the left action of (H and their 1 2 counterparts are the only elements that can differ betweenK(H 2 ) and (K(H)) 2 .
A. Diagonal elements : H + 1 H − 1 As mentioned previously, applying H − i to |(n 1 , n 2 ) does not give a symmetric superposition. In order to describe correctly the resulting state we need to introduce more detailed classes.
(S19) Using these new classes and identities, it is now possible to compute the matrix element forK(H + 1 H − 1 ) acting on |(n 1 , n 2 ) .
(S24) m 1 and m 2 can be viewed as random variables with a probability distribution proportional to #(n 1 , n 2 , m 1 , m 2 ), then Eq. (S24) is simply the expectation value E[m 1 m 2 ]. The equivalent element for (K(H)) 2 can itself be understood as a product of the expectation values where the relation between E[m i ] and the matrix elements ofK(H) is derived from Eqs. (S18-S19). The subspace energy variance for H + 1 H − 2 is then given by the difference between equations Eq. (S24) and Eq. (S25), which is equal to the covariance of m 1 and m 2 Up to a multiplicative factor, this quantity is equal to m1,m2 We will now make the assumption that N 2 −n 1 −1 = 0 and N 2 −n 2 −1 = 0 and demonstrate that this expression is identically zero. The cases n i = N 2 − 1 will be handled separately later.
For the purposes of this proof we will introduce an ordinary generating function for the elaborated class sizes, From Eq. (S27) we can derive a differential equation for this generating function, ∂ y1 ∂ y2 (w∂ w +x 1 ∂ x1 −1)(w∂ w +x 2 ∂ x2 −1) x 2 x 2 x 1 y 2 y 1 y 2 x 1 x 1 y 1 1 Figure S1. Finite state automaton diagram for the elaborated generating function. From this diagram we can read off the matrix elements of Eq. (S30).
where we have eliminated z by a change of variable in favour of w = z/(x 1 x 2 ). This equation is satisfied if and only if Eq. (S27) is equal to 0.
We will now derive a rational expression for G. Take a look at Fig. S1. Each state of this automaton is the state of the two sites on the (left) boundary of a partial configuration. Each arrow is an allowed gluing process. It takes the motif on the target of the arrow and glues it onto the left on the partial configuration producing a configuration with two more sites. The label on the arrows is a monomial factor which is picked up when following that gluing process. The x 1 and x 2 indeterminates mark when an excitation is added to the respective sublattices, whereas y 1 and y 2 are placed when we realise that we could have placed an excitation. The sum of all the execution paths of this automaton generates all the allowed configuration and marks them with the appropriate powers of the indeterminates. A process matrix for this automaton is the following, with which we can express the generating function as a Neumann series. The series clearly is convergent in some neighbourhood of z = x 1 = x 2 = 0 and y 1 = y 2 = 1 and so we may express it as where This differs by an additive constant from the definition given previously, but this detail is unimportant. It is a lengthy albeit straightforward calculation to verify that this indeed satisfies the differential equation Eq. (S29) as required. Now we return to the deferred cases of when either N 2 − n 1 − 1 = 0 or N 2 − n 2 − 1 = 0. The first term E[m 1 m 2 ] vanishes because in either of these cases, for one of the sublattices, it is not possible to insert an excitation, hence m 1 m 2 = 0. For the second term E[m 1 ]E[m 2 ] the same reasoning implies that one of #(n 1 + 1, n 2 ) and #(n 1 , n 2 +1) must vanish. Thus the proposition is shown inclusive of the previously excluded cases.

C. Total subspace energy variance
The only non-zero matrix elements of the subspace energy variance are the diagonal ones given by Eqs. (S21-S22). Taking their difference leads to for |(n 1 , n 2 ) . The total subspace variance can also be obtained by tracing over all basis states

SIII. Proof that M ⊆ K
We will work in the general setting where there is a collection of N sites, an assignment of them onto d subsets and a projection operation which is diagonal in the computational basis. The state space M are the projected 2-level coherent states and the vector space K is the span of the symmetric combinations of the computational basis states with fixed occupation on each of the subsets. This setting subsumes the cases of the coherentstate TDVP ansatz for the PXP model from Ref. 9 developed further in Ref. 11, and also the generalisations for Z d initial states and tree-tensor networks. 24 The states in M can be expressed using matrix product states with bond dimension D = 2. Similar to Ref. 11, we use the parametrisation where A • is for an excited site and A • for an unexcited site. The chain is split into d sublattices, each with two parameters θ j and φ j . This produces N/d identical unit cells, and the corresponding (unnormalised) quantum state is

S5
where the trace enforces periodic boundary conditions. This is equivalent to fixing sites j, j + d, j + 2d, . . . to cos θj 2 |• − ie iφj sin θj 2 |• , with the additional application of a projector that annihilates configurations with neighbouring excitations.
Up to this point, only the symmetric subspace with two sublattices (d = 2) has been considered. However, it is easy to see how to extend our previous results to d sublattices. The basis states |(n 1 , ..., n d ) N are then symmetric superpositions of product states with n i excitations on the i-th sublattice. As in Sec. SI, recurrence relations for the class sizes can be found using combinatorics arguments. Once theses are known, the matrix elements can be computed from them as it was done for d = 2. We will now demonstrate that, for any value of d, all states in M lie in the corresponding K subspace.
In Eq. (S38), for each configuration of σ's there is a product of N matrices of size 2 × 2. For the moment, suppose that there is at least one excitation in the configuration. Then this sequence can be segmented into blocks with the following pattern: the leftmost matrix corresponds to an excitation and all the others are vacancies. Because of the cyclic property of the trace, we can also assume that the leftmost matrix in Eq. (S38) corresponds to an excitation and rewrite the quantity inside the trace as where the indices of the angles are periodic in d, meaning that θ d+j = θ j . From Eq. (S37), it is easy to compute the product of two vacancy matrices Therefore, in each block the whole trail of vacancy matrices on the right can be reduced to a single vacancy matrix and a product of cosines. The product between an excitation matrix and a vacancy matrix can also be simplified as where B := 1 0 0 0 . Finally, it is easy to see from Eq. (S37) that the product of two excitation matrices gives 0. So this parametrisation guarantees that each block has at least one vacancy matrix on the right, and that Eq. (S39) can be written as This is a product of collinear (and hence commuting) matrices, and the effect of the trace on it is trivial. The last step is to count the number of blocks with a leftmost matrix with index j. For any configuration, it is simply given by the number of excitations in the corresponding sublattice. This means that each excitation in the sublattice j gives a factor of ie iφj sin(θ j /2). The other sites in this sublattice are free and give a factor of cos(θ j /2) each. This means that each configuration n 1 , . . . , n k has a prefactor of SIV. Proof that K = span M Section SIII establishes that M ⊆ K. We now claim that the relation between these two spaces is stronger, and that K = span M. For simplicity, in this section we will work with the complex parametrisation of M, (S46) This parametrisation, like the coherent-state parametrisation of section SIII, may fail to be injective. This happens for the PXP model when θ a = ±π/2 for some a. This is a defect with the space M, where M fails to be locally Euclidean, rather than with the parametrisations.
We will now be considering two ordered bases B and B where span B = K and span B ⊆ span M . The first basis B is defined as B n1,...,n d = #(n 1 . . . , n d ) N |(n 1 , . . . , n d ) N , (S47) These can be read to mean that each basis vector in B is a linear combination of those in B, and vice-versa, i.e span B = span B . Recalling that B is a complete basis for K and that every vector in B is contained in M, we find the containment, span M ⊇ span B = span B = K, which completes the proof.

SV. Bipartite entropy in the symmetric subspace
Computing the von Neumann bipartite entropy in the symmetric subspace is not straightforward due to the non-local definition of the basis states. Consider a spin chain of length N with periodic boundary conditions and cut it into two subsystems of size N L and N R , respectively. Let quantities related to the left part be denoted by the subscript L and the ones related to the right part by the subscript R. The simplest way to compute the bipartite entropy for this cut for a basis state |(n 1 , n 2 ) N would be to express it in the full Hilbert space and proceed from there. However this is suboptimal as it would strongly limit our ability to study large systems. Instead, the structure of the subspace can be exploited in order to write the states in the right and left subsystems in their respective symmetric subspaces. To simplify the computations, we only consider cases where the right and left subsystems are of even size. If nothing is specified, we consider the cut into two subsystems of equal size, meaning that the full system has size N = 4M, M ∈ N, and the two subsystems have size N L = N R = 2M .
While Let us denote by (n 1 , n 2 , a, b) N , a, b ∈ 0, 1 the class of states in (n 1 , n 2 ) N with the additional constraint that the leftmost site is defined by a and the rightmost site is by b. The convention used is that the leftmost site is unexcited if a = 0 and excited if a = 1, and the same hold for b and the rightmost site. Then |(n 1 , n 2 , a, b) N is the normalised symmetric superposition of all product sates in (n 1 , n 2 , a, b) N . Using these states, the relevant Schmidt decomposition for |(n 1 , n 2 ) N is where the sums on l i and r i run from 0 to N 2 , the sums on all a and b run from 0 to 1, and α depends on n 1 , n 2 and on all the summation indices. This decomposition includes terms that are forbidden by Rydberg constraints such as |( N L 2 , N L 2 , a L , b L ) N L . However we will see later that the coefficient α corresponding to these will naturally be zero.
To compute α, selection rules can be stated. The most obvious one is the conservation of the number of excitations on each sublattice, which tells us that l i + r i = n i . More formally, it means that α needs to include a term β n1,n2 l1,l2,r1,r2 = δ l1+r1,n1 δ l2+r2,n2 , where δ i,j is the Kronecker delta. The selection rule on the a and b is also simple: a L and b R cannot both be 1 and neither can b L and a R as this would correspond to neighbouring excitations. This can be expressed as The last step is to compute the contribution of each class (n 1 , n 2 , a, b) N . Each product state in (n 1 , n 2 ) N corresponds to a single state in (l 1 , l 2 , a L , b L ) N L ⊗ (r 1 , r 2 , a R , b R ) N R . So the only multiplicative factors that enter the equation are the normalisation ones. For each combination of the l i , r i , a, b the contribution to α is simply given by

S7
These class sizes can be computed using the same arguments as in section SI, keeping in mind that the two subsystems have open boundary conditions. To separate the contributions of the left and right parts we define γ N n1,n2,a,b = #(n 1 , n 2 , a, b) N .
This coefficient ensures that all states in the Schmidt decomposition respect the Rydberg constraint. Indeed, if a state is incompatible with it, its class size must be 0 and so is the corresponding γ. Using Eqs. (S52-S55), the coefficients α from Eq. (S51) can be written as From this Schmidt decomposition it is straightforward to compute the entanglement spectrum and the von Neumann entropy.

SVI. Revivals in large-N limit
The main feature of the PXP model is the revivals of the wavefunction that can be seen for some initial product states. This effect is the strongest when starting from the Néel state |Z 2 = |•••• . . . or the anti-Néel state |Z 2 = |•••• . . . . When starting in one of these states, after a time T the wavefunction comes back to itself. Furthermore, when starting in the state |Z 2 , after a half-period the wavefunction will be |Z 2 , and vice-versa. In the remainder of this section we will refer to this as state transfer. However, these revivals (and state transfers) are not perfect and decay with time. On the other hand, TDVP shows a periodic orbit going through both the Néel and anti-Néel states for an infinite system. In this section we study the finite-size dependence of revivals in the symmetric subspace.
To measure the quality of the revivals, the wavefunction fidelity and the entanglement entropy are used. Starting from the Néel state |( N 2 , 0) N = |Z 2 , the Schrödinger time evolution is generated by e −iHt . The fidelity of revival is then and the fidelity of state transfer is (S58) For the entanglement entropy we use the von Neumann bipartite entropy, as computed in Sec. SV. We also restrict our study to systems of size N = 4M, M ∈ N, and we cut them into two equally sized subsystems. From this we define the entropy as As both the Néel and anti-Néel states are product states in the full Hilbert space, their entanglement entropy is 0. From the exact PXP model, we expect the first state transfer to happen around t = 2.35 and the first revival around t = 4.70. 10 This can be seen in the entropy plot in Fig. S2. From this plot it is clear that the revivals get enhanced as N increases, but the quality of the first state transfer vary very little with system size. To obtain a more rigorous scaling, let us define the times of the first revival and state transfer, respectively, as This leads to their respective fidelity and entropy The scaling of these quantities with N can be seen in Fig. S3. The results for f 1 are consistent with the symmetric subspace having perfect revivals in the thermodynamic limit. However, the state transfer fidelity f 1/2 is clearly not converging towards 1.
To understand this behaviour we need to look at the different symmetry sectors as well as the different bands of states. With periodic boundary conditions, the relevant symmetry of the PXP model is translation. Because of the periodicity of the Néel state, it is only composed of eigenstates with momentum k = 0 or k = π. The scarred states of the top band belong alternatively to these two sectors, with the condition that the ground state always has k = 0. In the symmetric subspace, the quasimodes show the same behaviour with respect to these sectors.
Because our method allows us to access larger systems, more can be said about the scarred eigenstates.  Figure S3. Wavefunction fidelity and bipartite entanglement entropy at the first state transfer (f 1/2 , S 1/2 ) and revival (f1, S1). The fits shown are of the form y = a · 1/N + b. The data is for sizes up to N = 720 for f 1/2 , f1 and for sizes up to N = 560 for S 1/2 , S1.
For example, scaling the system size shows that there is a constant fraction of them with an equal energy spacing around E = 0. As N → ∞, the number of scarred quasimodes with a significant support on |Z 2 scales as O( √ N ). This means that, as we approach the thermodynamic limit, the only scarred eigenstates that have an non-zero overlap with the Néel state are evenly spaced.
Equal spacing of the energies with a difference ∆E implies perfect revivals with a period T = 2π/∆E. Without loss of generality, let us assume that the energy of the jth eigenstate is equal to j · ∆E. After half a period, the phases of the even and odd eigenstates are Because of the alternation of momentum, only the k = π states get a −1 factor as they are the odd ones. This exactly corresponds to the anti-Néel state as whereT is the translation operator. This means that if only the top band were relevant there would be perfect revivals and perfect state transfer in K. However, the overlap between |Z 2 and the non-scarred states is not negligible. In particular, we can explain the discrepancy between revivals and state transfers by looking at the second band. By "second band", we understand the band of states with the highest overlap with |Z 2 after the scarred states (see figure S4). The quasimodes in this band have a much lower overlap with the Néel state than the scarred ones, but they also display the alternation between the two momentum sectors. Near E = 0, the two top bands converge towards the same energy values (Fig. S5), but with the opposite momentum value. Hence, in this region there are two eigenstates for each energy value, one with k = 0 and the other with k = 1. Because of that, the previous argument for state transfer does not hold anymore. However this has no effect on the energy spacing itself, so perfect revivals are still possible.
To summarise, in the symmetric subspace revivals get better as system size increases. While the finite size scaling is consistent with revival becoming perfect in the thermodynamic limit, it appears that the state transfer is never perfect because of the presence of the second band. This is to be contrasted with TDVP where both revivals and state transfers are exact.

SVII. Path-integral construction
If we seek to model a particular quantum dynamics with a classical analogue, we need a good understanding of the relationship between the two. In the main text we stated that in a certain situation we can view the TDVP classical system as the dequantisation of the original dynamics. This is in the sense of path-integral quantisation and with the proviso that from the state space we can form a good resolution of the identity. The following subsection will provide an asymptotic bound on the fluctuations in S −1/2 µ . This bound establishes convergence between the transformed and untransformed classical systems in the thermodynamic limit, and shows how a measure can be constructed with which the untransformed system quantises naturally.

A. Path integral dequantisation
The quantum path integral comes in a number of conventional representations: the position or momentum representations, corresponding to the Schrödinger formulation, and coherent-state representation, corresponding to the phase space formulation. In each case, the quantisation procedure involves an integral over all admissible paths, including those with non-stationary action, which assigns to each path a weight given by a measure function and a phase provided by the classical action. Because our starting point is instead a quantum system, we will instead take the backwards view where the path integral is defined by a classical limit. We are also provided with a collection of states |x for which we assume that we can pick a measure µ to form a resolution of the identity, 1 =ˆdµ(x) |x x| .
The first term here gives another contribution to the path integral measure, and the second produces the kinetic term in the resulting Lagrangian, Therefore, in this scenario where we can assume Eq. (S69), the TDVP action describes a valid classical limit and this path-integral provides a way to quantise it to recover the full quantum dynamics. Many of the familiar manifolds of quantum states, such as spin coherent states and matrix product states come with a natural choice for this measure which is the Haar measure for a transitive group action. The constrained coherent states of M, however, are not of this type.

B. Fluctuation bounds for the frame operator
We begin by recalling the development of the transformed frame in the main text, before providing a proof for the asymptotic equivalence between the transformed and original frames. The first frame we considered involved the gauged parametrisation, Eq. (6), for M equipped with the measure dµ = N 2 + 1 4π 2 i=1,2 dθ i dφ i sin θ i cos 2 θ i 2 .
This measure was motivated with an intuitive picture of fuzzy distinguishability on M in the main text. From this the matrix elements of the frame operator S µ in the symmetric basis can be expressed as an integral. The integral over the azimuthal parameters φ a is oscillatory, ensuring that only the diagonal matrix elements are nonzero. The diagonal matrix elements are now computed, where n = n 1 + n 2 and the integration is performed by taking a substitution t i = cos 2 (θ i /2) to place it in the form of an Euler integral. For periodic boundary conditions, we arrive at the rational function which was quoted in the main text, where (x) (n) is the rising factorial or Pochhammer polynomial. 81 This is best thought of as a function of the densities, n a /N . Each linear factor in the numerator and denominator is a zero line or pole line, respectively, for the frame operator. For finite N and with the exception of the Néel states, these lines are outside of the domain of allowed densities. On the Néel states, the zero line (N/2 − n) and pole line (N/2 − n a ) intersect and cancel. In the vicinity of the Néel configurations, there are three 'corner' pole lines with fixed n a and three 'diagonal' zero lines with fixed n, all within a distance of order 1/N . As N increases the pole lines may become closer, but proportionally the zero lines also become closer. The scaling imparted by the zero lines to the residues of the approaching poles removes the divergence that would otherwise be caused by their approach. The additional cosine factors in the measure have the effect of adding two diagonal zero lines and one pole line passing-by each Néel state, and so without this there would be a residual pole at each Néel state, spoiling boundedness. Along the edge close to (N − n) = 0, there are no poles to cancel the zero lines so we can find matrix elements which are arbitrarily small as N increases. This provides a frame operator which is bounded and positive definite, even in the N → ∞ limit, except along the triangle edge. This last point is not too concerning because we know that the dynamics avoids this part of K, and is greatly preferable to leaving a simple pole on the most physically important states.
As emphasised in the main text, the frame operator found in this way is not suitable for quantisation. Its bounded and essentially positive-definite nature, however, leaves a measure which is qualitatively close to one in which a good resolution of the identity is found. Using the resulting frame operator we defined a transformed frame M, which produces a resolution of the identity from a congruence, This frame transformation potentially disturbs the states and obscures their physical meaning. We remedy this situation by approximating the effect of the frame transformation with a mere measure transformation, with the transformed measure, which we claim is asymptotically exact. This measure in the thermodynamic limit is shown in Fig. S6, where it is seen to be smooth away from the corners. This allows the classical system produced by the retraction onto the gauged constrained coherent states to quantise naturally to the quantum dynamics projected into K, as described in the previous subsection. The error in the approximation can be characterised by the quantity, This is the norm of the difference between the action of S −1/2 on a state |Ψ(θ, φ) and the action of multiplying by the state by the expectation value of that operator. It can be used to bound the difference between the measure-transformed frame operator S ν and the identity,