Emergent quantum state designs from individual many-body wavefunctions

Quantum chaos in many-body systems provides a bridge between statistical and quantum physics with strong predictive power. This framework is valuable for analyzing properties of complex quantum systems such as energy spectra and the dynamics of thermalization. While contemporary methods in quantum chaos often rely on random ensembles of quantum states and Hamiltonians, this is not reflective of most real-world systems. In this paper, we introduce a new perspective: across a wide range of examples, a single non-random quantum state is shown to encode universal and highly random quantum state ensembles. We characterize these ensembles using the notion of quantum state $k$-designs from quantum information theory and investigate their universality using a combination of analytic and numerical techniques. In particular, we establish that $k$-designs arise naturally from generic states as well as individual states associated with strongly interacting, time-independent Hamiltonian dynamics. Our results offer a new approach for studying quantum chaos and provide a practical method for sampling approximately uniformly random states; the latter has wide-ranging applications in quantum information science from tomography to benchmarking.

Analyzing the exact dynamics of general, strongly interacting quantum many-body systems is intractable using existing analytic and numerical tools. However, there is a widely used heuristic for understanding chaotic quantum dynamics: the eigenstates and eigenvalues of chaotic Hamiltonians have properties as if they were sampled from a random ensemble. This heuristic leads one to leverage statistical approaches, such as random matrix theory [1,2], to address many physical problems and has become a foundational principle in understanding chaos and thermalization in quantum systems [3][4][5][6]. Examples of this heuristic include Berry's conjecture [7] and the eigenstate thermalization hypothesis (ETH) [8][9][10], which have been supported by an overwhelming amount of numerical evidence [11][12][13].
While these statistical approaches hinge on the presence of random ensembles, realistic quantum systems are often described by a fixed Hamiltonian. Accordingly, energy eigenstates, eigenvalues, and states evolved by Hamiltonian dynamics are deterministic functions of a particular Hamiltonian without any randomness involved. Thus, it is of fundamental interest to establish a connection between isolated quantum systems and the emergence of random ensembles that dictate their statistical properties.
In this paper, we present a new perspective on the emergence of statistical behavior in chaotic quantum many-body systems: instead of imagining that a physical state is sampled from a random ensemble, we use a single many-body wavefunction to generate an ensemble of pure states on a subsystem; this ensemble goes beyond the reduced density matrix. Then we leverage statistical tools to analyze the degree of randomness of the local ensemble, and in turn characterize the global wavefunction.
Concretely, an ensemble of states can be generated from a single wavefunction by performing local measurements over only part of the total system. We consider a many-body system partitioned into a subsystem A and its complement B. Performing local measurements on B, we obtain exponentially many different pure states on A, each corresponding to a distinct measurement outcome on B. We call the set of pure states on A, along with the associated measurement probabilities, the projected ensemble on A (see Fig. 1a,b); see also [14]. To quantify the degree of randomness of an ensemble, we use a well-established notion from quantum information theory, namely quantum state k-designs [15,16]. An ensemble of pure states is said to form a quantum state kdesign if the probability distribution over the constituent states is indistinguishable up to k-th moments from the uniform distribution over the entire Hilbert space. Therefore, forming a higher k-design implies that the ensemble is more uniformly distributed over the Hilbert space.
Remarkably, using our new approach, we find that approximate k-designs arise from a variety of natural quantum many-body states. We establish two theorems showing that the projected ensemble coming from generic many-body quantum states forms an approximate k-design as long as the size of B is sufficiently larger than the size of A. Furthermore, a concurrent work [14] finds evidence that approximate k-designs emerge from projected ensembles in a Rydberg quantum simulator. We argue that this is a much more general phenomenon: we find strong numerical evidence across several models that approximate k-designs arise from both quantum Emergence of a universal quantum state ensemble from a single many-body wavefunction. a, A subsystem B of a pure many-body wavefunction |Ψ is measured in a fixed local basis; the remaining unmeasured qubits in A are in a pure state that depends on the measurement outcome on B. b, For quantum systems consisting of qubits, the measurement on B samples random outcomes, each characterized as a bitstring zB (binary numbers in the blue/yellow arrows). Different measurement outcomes zB occur with probability p(zB) and lead to distinct quantum states |ΨA(zB) , forming the projected ensemble E = {p(zB), |ΨA(zB) }. Right panel: the ensemble of pure states |ΨA(zB) (black arrows) are randomly distributed in the Hilbert space of A (black sphere), forming an approximate quantum state k-design. c, Examples of manybody wavefunctions whose projected ensembles form approximate quantum state designs include the typical output states of random unitary circuit evolution, quantum states obtained from quenched time evolution, and energy eigenstates of a chaotic Hamiltonian at infinite temperature. d, Illustration of minimal quantum state k-designs for a single qubit on the Bloch sphere. Forming higher k-designs requires a larger number of pair-wise non-orthogonal quantum states. In the limit of k → ∞, a k-design approaches the so-called Haar ensemble, which is the uniform distribution over all pure quantum states.
states obtained by quenched time evolution and energy eigenstates of chaotic Hamiltonians (see Fig. 1c). In the former case, we find that the degree of randomness in a projected ensemble continues to grow even after conventional local thermalization has occurred. Specifically, we establish numerically that states evolved for a longer amount of time form higher approximate k-designs. In the case of energy eigenstates, approximate k-designs emerge from states in the middle of the energy spectrum corresponding to effective infinite temperature. For finite-temperature eigenstates we observe that the projected ensembles converge to a universal ensemble that smoothly varies with respect to the energy density.
Our findings suggest that for a wide range of physically relevant many-body states, projected ensembles exhibit a universal form of randomness. This allows us to quantify the chaotic nature of Hamiltonian dynamics and to study the growth of complex, nonlocal correlations between a subsystem and its complement beyond the conventional paradigm of quantum thermalization [8-11, 13, 17, 18]. Post-thermalization physics in quantum many-body systems exhibits interesting quantitative and qualitative differences from its classical counterpart (e.g. [19]), and the current work clarifies and extends this paradigm. Furthermore, our work presents a simple protocol to efficiently produce an ensemble of random pure states from fixed, time-independent Hamiltonian dynamics that does not require highly fine-tuned controls.

I. QUANTUM STATE DESIGNS FROM PROJECTED ENSEMBLES
Studying properties of random quantum states is useful because they often encode universal phenomena found in nature. A standard approach is to consider the Haar ensemble, which is the uniform distribution over all pure states in a Hilbert space. While the Haar ensemble can be studied analytically using various statistical tools [20], it is extremely challenging to experimentally realize, as doing so requires an exponential amount of resources (such as the number of quantum gates or experimental operations) [21].
Instead, ensembles that form quantum state designs are considered because they mimic the Haar ensemble and can be efficiently realized in physical systems [22][23][24][25][26][27]. To build intuition for quantum state designs, let us consider quantum state designs for a single qubit (see Fig. 1d). For example, the ensemble of singlequbit states {|0 , |1 } with equal probabilities has the same mean (first moment) as the Haar ensemble: , where E Ψ∼E denotes averaging |Ψ over the ensemble E. We say that such a two-state ensemble forms a 1-design. However, the two-state ensemble does not form a 2-design because its second moment 1 2 (|0 0| ⊗ |0 0| + |1 1| ⊗ |1 1|) differs from that of the Haar ensemble. To form a 2-design for a single qubit, one can use four distinct non-orthogonal quantum states uniformly spread over the Bloch sphere (Fig. 1d). In general, for an ensemble to form a higher order design, it must be supported over a larger number of states [28].
Formally, given an ensemble E = {p i , |Ψ i } consisting of a set of wavefunctions |Ψ i and their probabilities p i , one can test whether the ensemble forms an approximate k-design by explicitly comparing its k-th moment against that of the Haar ensemble ρ Haar (which is defined by a continuous probability distribution over pure states). More precisely, we say an ensemble E is an ε- where · 1 denotes the trace norm. This definition means that the k-th moment of E is nearly indistinguishable from the k-th moment of the Haar ensemble up to a small error ε. It can be shown that an ε-approximate k-design is also an ε-approximate j-design for any j < k, and accordingly larger values of k indicate that an ensemble looks more uniformly random. Unlike the Haar ensemble, approximate designs arise in physical settings [22][23][24][25][26][27]. A canonical example is random unitary circuits [22,23,25], where a set of random two-qubit unitary gates in a certain geometric arrangement are sequentially applied to simple initial states (for an example, see Fig. 1c). Then the ensemble of resulting states over different choices of unitary gates forms an εapproximate k-design as long as the depth of the circuit is sufficiently large, scaling polynomially in k, log(1/ε), and the number of qubits N [22]. Similarly, there are a number of proposals to generate approximate designs based on time-dependent local Hamiltonian evolution [24,27].
In these examples, states are sampled from approximate k-designs by realizing many distinct, highly engineered quantum evolutions. By contrast, we show that approximate k-designs arise naturally from the projected ensemble of a single many-body wavefunction. More precisely, consider a many-body wavefunction |Ψ in the Hilbert space H for a bipartite system consisting of N A qubits in A and N B qubits in B. The projected ensemble for A is generated by performing projective measurements on all N B qubits in B in a local basis {|z B }, where a bitstring z B ∈ {0, 1} N B enumerates over all 2 N B measurement outcomes. After the measurement, with probability p(z B ) the system is described by the normalized and 1 A is the identity operator acting on qubits in A. This defines a projected ensemble: The ensemble consists of 2 N B states which are generally not pairwise orthogonal. The projected ensemble contains strictly more information than the reduced density matrix, since the density matrix is the first moment of the ensemble. We will call |Ψ the generator state of E Ψ,A . We note that Refs [29, 30] use a related construction to define localizable entanglement. We find that projected ensembles form approximate k-designs for many physically relevant many-body wavefunctions, including special, analytically soluble cases. Specifically, it can be shown that a generator state for an ε-approximate k-design can be efficiently prepared by a simple protocol, in which an initial product state is evolved by a Hamiltonian over a short time duration. The Hamiltonian can be time-independent and only has nearest-neighbor Ising interactions. This is possible because a universal resource state for measurement-based quantum computation can be prepared within a constant time independent of system size [31,32]. While this example is fine-tuned, we show that, in fact, generic many-body wavefunctions are good generator states for ε-approximate k-designs: Theorem 1. Let |Ψ be chosen uniformly at random from the Hilbert space H. The ensemble E A,Ψ forms an ε-approximate k-design with probability at least 1 − δ if Here Ω(·) denotes a lower bound up to a constant multiplicative factor and subleading terms. This theorem is proved in Appendix C, and establishes that all but a tiny fraction (of order ∼ 1/2 2 N B ) of the states in the Hilbert space are generator states for approximate k-designs if N B is asymptotically larger than k times N A . However, a quantum state randomly sampled from the entire Hilbert space is not so physical since such a state is extremely difficult to produce experimentally [21]. To this end, we also present the following theorem: Theorem 2. Let |Ψ be a state sampled from an ensemble on H that forms an ε -approximate k -design. Then the projected ensemble E A,Ψ forms an ε-approximate kdesign with probability at least 1 − δ if The proof is given in Appendix C and relies on higherorder concentration of measure results [33] as well as a polynomial approximation technique used in quantum algorithms for solving linear systems [34]. This theorem shows that if the generator state is complex enough, in the sense that it is a typical state from an ε -approximate k -design for small ε and large k [35], then the projected ensemble will well-approximate a quantum k-design. Theorem 2 has implications for ongoing experiments: if a single sample of an approximate design is experimentally realized, our theorem states that it can be used to generate ensembles forming approximate designs on Emergent quantum state designs from chaotic time evolution. a, Quenched dynamics under a timeindependent Hamiltonian H starting from an initial product state, |0 ⊗N , for a N -qubit system. The system is partitioned into two subsystems A and B with size NA and NB, respectively. At time t, a projective measurement in the local z-basis is performed on subsystem B, resulting in a specific outcome zB of length NB. b, Trace distances ∆ (k) between the k-th moments of the Haar ensemble and projected ensembles for an NA = 3 subsystem as a function of evolution time for various total system sizes N = 12, 14, . . . 24 (darker colors for increasing N ). Dashed lines are a phenomenological power-law fit, yielding a scaling of ∆ (k) ∼ t −1.2 . c, Design time τ k defined by the evolution time to achieve a trace distance of ∆ (k) = 0.02. We find that longer time evolution is required to form a higher approximate k-design. d, Late-time trace distances at t ≈ 10 3 . For a mixed-field chaotic Hamiltonian (circles), the late-time trace distances exhibit an exponential scaling with NB while they remain nearly constant for an integrable, non-ergodic Hamiltonian (crosses, time traces not shown). Dashed lines represent the trace distances ∆ its subsystems. Moreover, this protocol is hardwareefficient since it does not require fine-tuned controls to produce many different quantum states, and could lead to various useful applications in quantum information science [22,[36][37][38][39][40].
At a conceptual level, Theorem 2 establishes that a large class of states which can be efficiently prepared are good generators of k-designs. This raises the possibility that, even in natural chaotic quantum systems, approximate k-designs may arise from projected ensembles.

II. QUANTUM STATE DESIGNS FROM CHAOTIC SYSTEMS
Motivated by the above results, we numerically investigate projected ensembles that arise from chaotic Hamiltonian dynamics. We start with the paradigmatic example of the 1D quantum Ising spin system with mixed fields (QIMF), described by the Hamiltonian where N is the number of spins, σ µ j with µ = x, y, z are the Pauli operators for a spin at site j, J is the strength of Ising interactions, and h x and h y are the strengths of the longitudinal and transverse fields, respectively. In the absence of the longitudinal field (h x = 0), the Hamiltonian can be mapped to an integrable model of non-interacting fermions via the Jordan-Wigner transformation, leading to non-ergodic dynamics. However, for any non-zero longitudinal field (h x = 0), the Hamiltonian is ergodic, and its eigenvalues and eigenvectors are expected to have properties consistent with ETH predictions. This has been explicitly checked for a specific parameter set (h x , h y , J) = (0.8090, 0.9045, 1) [41] which we adopt for our study. We note that we do not find a qualitative differences in our results when using nearby parameter values.
We first consider the many-body state |Ψ(t) = e −iHQIMFt |Ψ 0 resulting from time evolution of the initial state |Ψ 0 = |0 ⊗N (Fig. 2a). Here |0 j and |1 j are the eigenstates of σ z j with eigenvalues +1 and −1, respectively. The initial product state, |Ψ 0 , has zero energy expectation value with respect to H QIMF , corresponding to the total energy of an infinite temperature state. This implies that local subsystems will relax to an infinite-temperature ensemble after a local thermalization time [8][9][10].
At any time t, the projected ensemble for a subsystem A is obtained by simulating projective measurements on the rest of the N B = N − N A qubits in the local zbasis (Fig. 2a). 1 In order to check if the projected ensemble forms an approximate k-design, we compare the k-th moment of the ensemble, ρ (k) E in Eq. (1), to that of the uniform ensemble ρ (k) Haar using the trace distance . For a fixed subsystem size of N A = 3, we numerically compute the trace distance ∆ (k) up to k = 4 as a function of time for various N B (Fig. 2b). In all cases, ∆ (k) decays in time, following a phenomenological power-law scaling ∆ (k) ∼ 1/t 1.2 , until it saturates to a value that is governed by finite size effects. We note that the extracted scaling exponent is non-universal and Hamiltonian-dependent. The saturation value of ∆ (k) decreases exponentially with N B (circles, Fig. 2d), exhibiting the scaling ∆ (k) ∼ 1/ √ 2 N B . If the power-law scaling persists at larger system sizes then ∆ (k) may decrease over an exponentially long time scale.
As a comparison, we also present ∆ (k) em which is the trace distance between the k-th moments of the Haar ensemble and the empirical Haar ensemble consisting of 2 N B states sampled uniformly at random on A (dashed lines, Fig. 2d). Since the empirical ensemble asymptotically approaches the Haar ensemble in the limit of infinite samples, ∆ (k) em is determined only by statistical fluctuations associated with having a finite number of quantum states. Remarkably, we find that the projected ensemble obtained from the quench dynamics shows a trace distance almost identical to that of the empirical ensemble of the same size, suggesting that the former is as uniformly random as the latter. By contrast, repeating similar calculations for the integrable model (h x = 0), we observe qualitatively different behavior where the trace distance to the Haar ensemble is much larger than in the non-integrable case (crosses, Fig. 2d). Furthermore, there is no appreciable dependence on system size. This is expected, since integrable systems do not locally thermalize and instead relax to a generalized Gibbs ensemble, and hence will not form 1-designs at effective infinite temperature [42].
Given the emergence of k-designs in asymptotic regimes for a chaotic Hamiltonian, it is natural to ask how long it takes for a subsystem to achieve an approximate design up to a small, fixed precision. To this end, we introduce a design time τ k defined as the time at which ∆ (k) becomes smaller than a certain fixed threshold ε. For a chosen threshold ε = 0.02, we find that the formation of higher k-designs requires longer time evolution (Fig. 2c). This observation is consistent with the idea that typical quantum states from higher k-designs are more complex and hence more difficult to prepare [35, 43, 44].
Next we investigate the properties of energy eigenstates. We repeat a similar analysis as above by replacing the time-evolved state |Ψ(t) with an energy eigenstate |E i of H QIMF . Figure 3 shows the trace distance ∆ (k) as a function of energy E i for a projected ensemble generated from |E i . We find a sharp dip at zero energy density corresponding to infinite temperature, which signifies the emergence of approximate quantum state designs ( Fig. 3a). At zero energy, ∆ (k) decreases exponentially as a function of system size for k = 1, 2, 3. However, the decay rate is slightly slower than in the case of quenched dynamics (Fig. 3b).
In addition to the QIMF, we have also studied two other ergodic Hamiltonian models in order to corroborate the universality of our findings (see Appendix B). Specifically, we consider a system of random, all-to-all coupled spin-1/2 particles as well as hard-core bosons with random all-to-all hoppings with particle number conservation. In both cases, we observe excellent convergence of projected ensembles to approximate k-designs. In the latter case, the measurement outcomes in B and corresponding pure quantum states A are strongly correlated owing to the particle number conservation symmetry; hence a naïve approach based on Eq. (5) does not lead to approximate k-designs. We instead introduce symmetryresolved projected ensembles by grouping certain subsets of measurement outcomes from B (see Appendix B for details); this does lead to approximate k-designs.
For chaotic Hamiltonians, projected ensembles forming k = 1 designs can be anticipated from the standard picture of quantum thermalization since the first moment of a projected ensemble simply corresponds to the reduced density matrix of a subsystem. The reduced density matrix approaching the first moment of the Haar ensemble, i.e. the maximally mixed state, follows from local thermalization at infinite temperature. However, the convergence of higher moments k ≥ 2 of the projected ensemble to higher k-designs is nontrivial and surprising. Such Universal ensemble at finite temperature. a, Trace distances between the second moments of projected ensembles generated from a pair of energy eigenstates at Ei and Ej for NA = 3 subsystems. We plot the pairwise distances ∆ (2) ij for every pair of eigenstates |Ei , |Ej , computed for system size N = 11. The distances are minimized when Ei ≈ Ej. The plot suggests the existence of a universal ensemble that depends smoothly on the energy density. The black dashed line indicates a cut defined by (Ei + Ej)/N = −0.6. b, Trace distances plotted as a function of energy density difference, |Ei − Ej|/N , along the black dashed line in a for various systems sizes N = 9, 10, . . . , 15 (darker colors for increasing N ). Inset: the distances at zero energy difference Ei = Ej decay exponentially with system size. The trace distances from the projected ensemble are comparable to those from the finite-size ensemble, ∆ (2) em , of the empirical Haar distribution (dashed line). Such an exponential scaling suggests the existence of a universal random ensemble at finite temperatures.
convergence cannot be explained by ETH alone, and suggests a new form of emergent randomness beyond the conventional framework of quantum thermalization.
A natural next step is to generalize quantum state designs (or Haar ensembles) to a finite-temperature setting. However, we are unaware of any appropriate analogue of the Haar ensemble at finite temperature. Such an ensemble, if it exists, would generally depend on the system Hamiltonian, and its first moment should (approximately) be a thermal state. While explicitly identifying properties of such an ensemble is an interesting future direction, here we find numerical evidence that such an ensemble exists. In Figure 4, we compute projected ensembles for all energy eigenstates of H QIMF and present the pairwise distances where ρ (2) i denotes the second moment of a projected ensemble generated from an eigenstate |E i . We find that ∆ (2) ij is a smooth function of energy up to small fluctuations, suggesting that the projected ensembles smoothly vary with energy as well (Fig. 4a). Also, ∆ (2) ij is minimized when the energy difference |E i − E j | is small (Fig. 4b), and in this regime ∆ (2) ij decreases exponentially with system size (Fig. 4b inset). These observations suggest that the second moment of the projected ensemble is indeed universal even at finite temperatures.

III. DISCUSSION AND OUTLOOK
The random quantum state ensembles considered in this paper are qualitatively different from conventional ones in quantum statistical mechanics, such as the microcanonical and canonical ensembles. These latter ensembles are fully specified by their corresponding density matrices (their first moment) and are used to evaluate expectation values of observables.
Our formalism concerns more general statistical properties (such as higher moments) of an ensemble with a large number of pure states that are generally pairwise non-orthogonal. As such, a projected ensemble encodes additional information about a subsystem. For example, one can ask the following information-theoretic question: how much information (i.e., classical bits) is required in order to specify the full wavefunction of a random sample from a projected ensemble? If the projected ensemble were to be the Haar ensemble, an exponential number of bits would be required. By contrast, for an ensemble which is uniformly distributed only over computational basis states of A, a linear number of bits suffice to specify a sample. Notice that the former and latter ensembles produce the same density matrix, namely the maximally mixed state. Therefore, projected ensembles provide a novel framework to analyze the information content associated with quantum states of a subsystem and their relation with the remainder of the system.
Our findings open up a number of new directions in quantum chaos, thermalization, and quantum information. In particular, our numerical results demonstrate that for an initial product state evolved by a chaotic Hamiltonian, the largest k-design attained by the projected ensemble grows as a function of time; moreover, this growth persists significantly past the local thermalization timescale. While it is not presently clear how long this growth will persist, the observed growth is clearly diagnostic of the sustained development of nonlocal correlations after local thermalization has occurred (e.g. [19,45]). There may even be connections to the quantum complexity of a state evolving by chaotic dynamics, which likewise grows long after thermalization has occurred [35,46]. It would also be interesting to fully generalize the above to projected ensembles at finite temperature, in the presence of symmetries, (quasi-)integrability, or strong disorder resulting in localization [11][12][13].
In quantum information science, quantum states designs are valuable resources in many applications. Our work establishes projected ensembles as a practical way of sampling states from approximate designs using natural Hamiltonian evolutions of existing quantum simulators without fine control. Further, our work could lead to novel experimental quantum tomography protocols [40], cryptographic protocols for hiding information [22], the design of unforgeable quantum encryption [36], and also new methods for quantum device verification [39,47]. Indeed, a parallel work [14] uses projected ensembles to devise and implement a novel benchmarking protocol.
Finally, an important question at the intersection of computer science and quantum many-body physics is whether the computational complexity of simulating natural chaotic dynamics is beyond the capability of contemporary classical computers. In other words, can quantum supremacy tests be performed using a fixed chaotic Hamiltonian with analog quantum simulators? Existing sampling-based quantum supremacy protocols heavily rely on certain statistical and computational properties of state ensembles formed by applying random unitaries to a fixed state [48-50]. The projected ensemble emerging from generic quantum dynamics may also possess the requisite properties and our work could lead to a new approach to quantum supremacy using analog quan-tum simulators. AFOSR YIP (FA9550-19-1-0044 To quantify the degree of randomness of our projected ensembles, we computed the trace distance ∆ (k) between the k-th moments of the projected ensemble and the Haar ensemble: Here ρ (k) Haar is the k-th moment averaged over the Haar ensemble. For a Hilbert space H with dimension d, it has the form [20]: .
Here, S k is the symmetric group on k elements, and π ∈ S k is an element of the group. Perm H ⊗k is a representation of S k on H ⊗k which permutes the tensor factors according to The inverses in the subscripts are chosen so that Perm H ⊗k (π)·Perm H ⊗k (π ) = Perm H ⊗k (π•π ) (i.e., the representation is a homomorphism of S k , and not an antihomomorphism). It is readily checked that Eq. (A2) can be written as where Π k is simply the projector onto the symmetric subspace of H ⊗k (that is, the invariant subspace under S k ); this subspace has dimension d+k−1 k .

Finite sampling error from the empirical Haar ensemble
In the main text, we constructed projected ensembles of size 2 N B , with N B the size of the complement subsystem. In Figs. 2d, 3b, and 4b, we compared the system size scaling of ∆ (k) against the trace distance ∆ (k) em of the empirical Haar ensemble, an ensemble formed by sampling from the Haar ensemble 2 N B times. In particular, we have This comparison is made to estimate the degree to which the error ∆ (k) in our projected ensemble is due to its finite size. We estimate ∆ i σ z j terms so that our initial state |Ψ 0 = |0 ⊗N has zero energy, and accordingly is regarded to be at infinite temperature.
As shown in Fig. 5a, we see excellent convergence to k-designs for our projected ensembles constructed from time-evolved quenched states as well as eigenstates. The late-time ∆ (k) values are very close to ∆ (k) em . The average ∆ (k) values for infinite-temperature eigenstates are also close to ∆ (k) em , and show clear exponential decay (Fig. 5b). Notably, we do not average over disorder realizations: each data point and time series are computed with fixed disorder realizations. This provides additional support for our conjecture that ergodic Hamiltonian systems generate approximate k-designs via time-evolved states and eigenstates at infinite temperature.

Random hopping model with U(1) symmetry
Having established that chaotic models such as the QIMF and the all-to-all random coupling model provide projected ensembles which converge to k-designs, we next ask whether quantum models with symmetries do so as well. We specifically study a random hopping model where J ± ij are random variables drawn from i.i.d. normal distributions: J ± ij ∼ N (0, 1/N ). Like the all-to-all random coupling model, this model is also expected to exhibit chaos. However, the specific interaction terms are chosen such that this model has a U (1) symmetry: there is conservation of the total magnetization S z tot. = j S z j . The model can equivalently be viewed as describing hard-core bosons with random complex all-to-all hopping amplitudes.
Any bitstring in the z-basis is a suitable infinite-temperature quench state. For N even, we choose the initial state |Ψ 0 = |0101 · · · 01 which lies in the S z tot. = 0 sector. Under time evolution by our chaotic Hamiltonian, the wavefunction is ergodic in the S z tot. = 0 Hilbert space. Naïvely forming the projected ensemble, we do not find convergence to an approximate k-design. This is shown in the light curves in Fig. 5c. With increasing system size, ∆ (k) for k = 1, 2, 3 saturate at a non-zero value. This is because of large correlations between the bitstrings z B and the projected state |Ψ A (z B ) : if z B has total magnetization s B , |Ψ A (z B ) necessarily has total magnetization S z tot. − s B due to the global conservation law. The Hilbert space H A naturally decomposes into a direct sum of multiple sectors enumerated by the magnetization s A . Accordingly, the projected ensemble now produces multiple ensembles, one for each sector.
In order to properly account for the U (1) conservation law, instead of projecting onto all strings z B , we post-select a subset of all strings: bitstrings with fixed total magnetization, e.g. s B = 1/2. The projected states |Ψ A (z B ) will also have fixed magnetization, and we then ask whether this ensemble forms a k-design in the subspace of H A with magnetization s A = S z tot. − s B . In our numerical examples, we present results with N A = 5 qubits and s A = −1/2, s B = 1/2. The relevant subspace of H A has dimension 10, far smaller than 2 N A = 32. Fig. 5c,d show the results of our symmetry resolution. We now see excellent convergence towards a k-design. As with the random coupling model, the late-time ∆ (k) values are very close to ∆ (k) em . In order to make a fair comparison, we plot the late-time ∆ (k) against an "effective N B ", defined as log 2 of the number of post-selected strings z B . ∆ (k) em is also computed by sampling the same number of times.
Using this post-selection procedure, we find excellent k-design convergence for projected ensembles from eigenstates near E = 0. Our results indicate that our basic approach remains valid for chaotic models with additional symmetries if these symmetries are properly addressed.
Appendix C: Proof of main theorems

Proof sketch for the main theorems
We begin with a sketch of the proof of Theorem 1. We will use streamlined notation for clarity. Recall that for a generator state |Φ on H, the projected ensemble It will be convenient to define the unnormalized state We note that in our new notation, we can write The key to establishing Theorem 1 is understanding the random tensor A(|Φ ) when |Φ is a Haar-random state.
Here, the Haar-random state is a normalized vector chosen uniformly at random in 2. The concentration of A(|Φ ) around its expectation. This will tell us that with high probability, A(|Φ ) is close to its expectation.
Let us discuss these two ingredients in more detail. Using the fact that the uniform measure (i.e., the Haar measure) on the complex sphere is invariant under any unitary rotations, we can show that p z and |Φ z Φ z | are independent random variables. Furthermore, |Φ z Φ z | is a uniform random vector from a d A -dimensional complex sphere. Hence, we have which is the k-th moment of the uniform ensemble. This means that the expectation of the k-th moment of the projected ensemble of a randomly selected many-body wavefunction |Φ reproduces the k-th moment of the Haar ensemble in d A complex dimensions. Hence the expectation E Φ∼Haar(d) A(|Φ ) is exactly equal to the desired quantity.
To understand the concentration of A(|Φ ) around its expectation, we recall the well-known result that the uniform distribution over a high-dimensional sphere has a very sharp concentration. For any Lipschitz function on the sphere, the fluctuation around its expectation value is small and the probability of having a large fluctuation decays exponentially. This is known as Levy's lemma and allows us to upper bound the probability that A(|Φ ) is far from its expectation. Together with the expectation identity, we have Prob Φ∼Haar(d) (C8) The asymptotic relation in Theorem 1 follows immediately from the above probabilistic statement. The detailed proof of Theorem 1 is given in Appendix C 3. The proof of Theorem 2 is very different from Theorem 1. First of all, the expectation of A(|Φ ) in Theorem 1 is computed using the invariance property of the Haar measure, which does not hold for a measure that only forms a state design. Furthermore, Levy's lemma only holds for the Haar measure, so we cannot resort to Levy's lemma to control its statistical fluctuations. To prove Theorem 2, we make use of a technique used in the context of solving linear systems on a quantum computer [34]. The key idea is to add a modulating function that approximates A(|Φ ) by a polynomial function in |Φ . In particular, where the modulating function is given by Using the binomial expansion, we can check that B(|Φ ) is a polynomial function. Furthermore, µ k,b (s) is very close to one when s = d B Φ z | Φ z is around one. Taking b larger allows us to better approximate the constant function µ k,b (s) = 1, which corresponds to the target expression A(|Φ ). A visualization of µ k,b (s) can be found in Fig. 6a. If |Φ is sampled from the uniform measure on the quantum state space, then s = d B Φ z | Φ z will concentrate around one; see Fig. 6b. So for |Φ sampled from the uniform measure, we can show that A(|Φ ) ≈ B(|Φ ). However, |Φ is sampled from a quantum state design, so there is no guarantee that s = d B Φ z | Φ z should concentrate around one.
To address this, we utilize the following result proved in Lemma 5 where R(|Φ ) is a polynomial function given by Hence, the approximation error can be upper bounded by a polynomial function that we have better control over. After introducing the key quantities, the proof proceeds by bounding the error between (i) A(|Φ ) where |Φ is sampled from a state design, and (ii) the expectation of A(|Ψ ) where |Ψ is from the Haar measure. We use B(|Φ ) as an intermediate point of comparison: 6. Visualization of the modulating function and the concentration of dB Φz| Φz . a, The modulating function µ k,b (s) defined in Eq. (C11) for varying values of b from 2 to 16. Darker colors correspond to higher b. We fix k = 2 in the plot. b, The concentration of s = dB Φz| Φz for the Hilbert space dimension of subsystem A, i.e., dA = 2 N A = 16, when |Φ is sampled uniformly at random from the complex sphere. We sample 10000 different vectors |Φ to generate the histogram. The black solid line is a kernel density estimation fit of the histogram. If dA is larger, then s is more concentrated.
The first term can be upper bounded by R(|Φ ) and we can use the fact that state designs approximate the expectation of any polynomial function under the Haar measure. Accordingly, In the second term, we can apply a similar idea by upper bounding the 1-norm with the 2-norm and utilizing the fact that B(|Φ ) − E Haar A(|Ψ ) is a polynomial function in |Φ . This gives We then approximate B(|Φ ) by A(|Φ ) in the above expression, where the error can be upper bounded by R(|Φ ) 2 . This gives the approximate relation Using these steps, all our expressions contain only expectations over the Haar measure. More precisely, The first two terms are small because 1. µ k,b (s) is close to one when s is close to one.
2. If |Φ is sampled from the Haar measure, then s = d B Φ z | Φ z will concentrate around one.
See Figure 6 for visualizations. The third term is small due to Theorem 1. Finally, we can resort to Markov's inequality to show that with a probability of at least 0.99, we have that A(|Φ ) − E Ψ∼Haar A(|Ψ ) 1 is small. The full proof of Theorem 2 is given in Appendix C 4.

Useful identities
We collect here several useful identities that we will leverage throughout the remainder of the Appendix. The first identity is Taking the trace, we find The related identity will likewise be useful.

Proof of Theorem 1
The main ingredient to prove Theorem 1 is Levy's lemma. Levy's lemma states that if a function is Lipschitz continuous, then it will concentrate around it's expectation value if the variable is sampled uniformly from a highdimensional complex sphere.
Before applying Levy's lemma, we will need to define a function of interest. In particular, we will consider individual entries in the main quantity z∈{0,1} N B p z |Φ z Φ z | ⊗k . Let {|i } i∈{0,1} tN A be the standard basis on H ⊗k A , and write |i = t k=1 |i (k) where each |i (k) ∈ H A . Similarly write |j = t k=1 |j (k) . Consider the function f ij : S 2d−1 → R defined by A nice property of this function is that it is Lipschitz continuous and hence Levy's lemma can be applied.

Lemma 2 (Lipschitz constant).
We have Proof. Since our f ij is differentiable, we can choose any η such that We have where we have simply explicitly evaluated the derivatives and used the triangle inequality for the 2-norm. Let us bound each of the term terms on the right-hand side in turn. For the first one, we have Eq. (C32) is less than or equal to Since M † p M ∞ ≤ 4, the above is ≤ 2t. Now we bound the term in the last line of Eq. (C30). Letting , the term can be written as Using the bound we can upper bound Eq. (C34) by Putting our bounds together, we have which is the desired result.
While the above bound on the Lipschitz constant, along with Levy's lemma, guarantee that each function f ij concentrates around the its expectation value, the expectation value of f ij has not yet been computed. The following lemma shows that the expectation value of f ij is closely connected to the quantum state k-design on subsystem A.
Lemma 3 (Expectation value identity). We have the identity Therefore the expectation of f ij is given by Furthermore, from a technical lemma given below, we have that |Φ z is independent from Φ z | Φ z . Using this fact and the linearity of expectation, we compute as desired.
During the evaluation of the expectation value, we used a property that the normalized state |Φ z Φ z | and the corresponding probability p z are independent random variables. This is proven in the following lemma.

Lemma 4.
If |Φ is a Haar-random state on H, then the random variables |Φ z Φ z | and p z are independent.
Proof. Define the map P z (|Φ ) := Φ|P z |Φ , and let R z be the map taking |Φ to the normalized state |Φ z , Let U A be a Haar-random unitary operator on H A and let U = U A ⊗1 B . Because U is unitary and |Φ is Haar-random, the state U |Φ is also Haar-random. Also note that where the equivalence is in the sense of random variables. Therefore, for any functions F and G, To pass from the second line to the third, we used the fact that for any state |Ψ on H A , the state U A |Ψ is an independent Haar-random state on H A . This shows that R z (|Φ ) and P z (|Φ ) are independent, as desired.
With the above lemmas, we now proceed with the proof of Theorem 1.
Proof of Theorem 1. Let us consider the function f ij (|Φ ) defined in Eq. (C26). Because we have bound the Lipschitz constant η ≤ 4t − 2, we can leverage Levy's lemma given in Lemma 1 to obtain Performing a union bound and rescaling ε → ε/d 2k A , we have (C45) By comparing to the definition of f ij in Eq. (C26) and using Lemma 3 to obtain the expectation value of f ij , the above is equivalent to the following concentration inequality where A entrywise, 1 = i,j |A ij |. Finally, using A entrywise, 1 ≥ A 1 and applying Lemma 3, we find We can see that if we have then the ensemble E Φ,A forms an ε-approximate t-design with probability at least 1 − δ. Recall that d A = 2 N A , d B = 2 N B , hence taking a logarithm on both side of the above condition gives rise to which is the result stated in the main text.

Proof of Theorem 2
For convenience we restate Theorem 2 in the following.
Theorem 3 (Restatement of Theorem 2). Let |Ψ be a state chosen from an ensemble on H that forms an εapproximate k -design. Then the projected ensemble E A,Ψ forms an ε-approximate k-design with probability at least In the following subsections, we will begin with a discussion of a generalized Levy's lemma that provides sharper concentration for quadratic functions. We will then give a general structure of the proof and provide the detailed proof of several technical lemmas afterwards.

a. Higher order concentration
In the proof of Theorem 1, Levy's lemma plays a crucial role in establishing the desired statement. Levy's lemma tells us about a sharp concentration when the random variables are sampled from a high-dimensional sphere. We will continue to utilize concentration on high-dimensional spheres for the proof of Theorem 2. However, the original statement of Levy's lemma does not provide good bounds for Theorem 2. We will instead make use of a higher order variant of Levy's lemma. In this section we recall a higher order variant of Levy's lemma that was established in [33], and use it to provide a simple proof of a concentration inequality for the quantity Φ z | Φ z . This higher order concentration results provided in [33] is a useful tool for obtaining concentration inequalities for polynomial functions on measures satisfying a log-Sobolev inequality, and may prove useful in other problems arising in quantum information theory. For this work we only need the following result from [33]: and Hess f (θ) ∞ ≤ 1 for all θ ∈ S N −1 then As a corollary, we obtain the following result for quadratic forms.
Corollary 1 (Concentration for quadratic forms). Let Q : R N → R N be a linear map satisfying Q 2 ≤ 1 2 and Q ∞ ≤ 1 2 , and let Proof. Let f (x) = x T Qx − a 0 . Then ∇f (x) = 2Qx, so Let {u j } N j=1 be an orthonormal basis of eigenvectors of Q T Q with eigenvalues λ 2 j . Then the latter integral is equal to If Q 2 ≤ 1/2 then f satisfies the condition (C54). Moreover Hess f = 2Q, so Hess f ∞ = 2 M ∞ . We can therefore apply Theorem 4 to f to obtain (C56).
We only need to use Corollary (1) in the case that Q is an orthogonal projection.
Corollary 2 (Concentration for quadratic forms with projectors). Let V ⊂ R N be a subspace of R N with dimension m, and let P V be the orthogonal projection onto V . Then and in particular, for θ sampled uniformly from the sphere, Proof. Observe that P V 2 2 = m and P V ∞ = 1, so the operator (2m) −1/2 P V satisfies the conditions of Corollary 1, from which (C59) follows. Then (C60) follows from Markov's inequality applied to (C59).
Specializing to the case that |Φ is a Haar-random quantum state on H and P V is the projector P z = 1 A ⊗ |x x|, we deduce an exponential concentration inequality for Φ z | Φ z .

Corollary 3 (Concentration for probability in a projected ensemble). We have the bound
The proof is based on the idea of modulating the target expression, which is a rational function, with a high-degree polynomial to form a polynomial function. More precisely, we will be interested in the following two functions: where r, b > 0 are parameters for tuning the approximation of polynomial function B(|Φ ) to the target expression A(|Φ ), which is a rational function. Because the 1 is cancelled in the binomial expansion it is not hard to see that B(|Φ ) is indeed a polynomial function in the real and imaginary parts of |Φ . We will consider |Φ to be sampled from an (ε , k )-design. We will use the basic Markov inequality to bound the concentration. A higher order concentration inequality can also be used but would require k to be higher. The central quantity in Markov inequality is the expectation value of the error (C65) We will use B(|Φ ) as a surrogate to obtain an upper bound on the above quantity. This is because B(|Φ ) is a polynomial function rather than a rational function in the real and imaginary parts of |Φ . A triangle inequality gives the following bound We can now analyze the two terms independently by using properties of quantum designs. For the first term, we will prove in Lemma 5 that the following inequality holds.
Therefore we will define a polynomial function which is an upper bound on the approximation error between A(|Φ ) and B(|Φ ). Because a quantum design approximates any polynomial function, we have the following upper bound where the exact expression of E 1 is given in Lemma 6. We will also apply a similar philosophy for bounding the second term in Eq. (C67) by first upper bounding the term by a polynomial function that turns expectations over designs into expectations over the Haar measure: The first line follows from the relation between 1-norm and 2-norm. The second line applies Jensen's inequality. The third line follows from the observation that d k A B(|Φ ) − E Ψ∼Haar A(|Ψ ) 2 2 is a polynomial function in |Φ . Therefore, we can turn the expectation over a quantum state design to one over the Haar measure, which incurs a small error of ε ×E 2 . The exact expression of E 2 is given in Lemma 7. We will now upper bound E Haar B(|Φ ) − E Ψ∼Haar A(|Ψ )   (C77) The first line follows from the triangle inequality and (a + b) 2 ≤ 2(a 2 + b 2 ), ∀a, b ∈ R. The second line follows from the relation between the 1-norm and 2-norm. The third line uses the inequality given in Eq. (C68), which is proved in Lemma 5. We can now combine Eq.'s (C67), (C71), (C74), and (C77) to find  . Along with Lemma 6 that bounds ε E 1 , Lemma 7 that bounds ε E 2 , we can show that for any ε > 0 as long as we choose r = 2 and b to be an even integer with b = Ω (N B + log(1/ε)) , (C80) N A = Ω (log(N B ) + log(k) + log log(1/ε)) , (C81) k = Ω (bk) , (C82) log(1/ε ) = Ω (log(1/ε) + k(bN B + log k + N A )) , (C83) we can obtain the following upper bounds Together with Markov's inequality, we have the following concentration result: for any ε,ε > 0. Using Lemma 3 established in the proof of Theorem 1, we have  Proof. Recall the following definitions: We note the definition of trace norm · 1 for Hermitian matrices: Hence, we have This concludes the proof.
where the error term is given by In particular, if we choose r = 2 ≤ d B and recall that d B = 2 N B , then for any ε > 0 as long as k = Ω (bk) , (C105) log(1/ε ) = Ω (log(1/ε) + bkN B ) , we have the following upper bound on the error term: ε E 1 ≤ ε.