Preparing quantum many-body scar states on quantum computers

Quantum many-body scar states are highly excited eigenstates of many-body systems that exhibit atypical entanglement and correlation properties relative to typical eigenstates at the same energy density. Scar states also give rise to inﬁnitely long-lived coherent dynamics when the system is prepared in a special initial state having ﬁnite overlap with them. Many models with exact scar states have been constructed, but the fate of scarred eigenstates and dynamics when these models are perturbed is diﬃcult to study with classical computational techniques. In this work, we propose state preparation protocols that enable the use of quantum computers to study this question. We present protocols both for individual scar states in a particular model, as well as superpositions of them that give rise to coherent dynamics. For superpositions of scar states, we present both a system-size-linear depth unitary and a ﬁnite-depth nonunitary state preparation protocol, the latter of which uses measurement and postselection to reduce the circuit depth. For individual scarred eigenstates, we formulate an exact state preparation approach based on matrix product states that yields quasipolynomial-depth circuits, as well as a variational approach with a polynomial-depth ansatz circuit. We also provide proof of principle state-preparation demonstrations on superconducting quantum hardware.

Quantum many-body scar states are highly excited eigenstates of many-body systems that exhibit atypical entanglement and correlation properties relative to typical eigenstates at the same energy density. Scar states also give rise to infinitely long-lived coherent dynamics when the system is prepared in a special initial state having finite overlap with them. Many models with exact scar states have been constructed, but the fate of scarred eigenstates and dynamics when these models are perturbed is difficult to study with classical computational techniques. In this work, we propose state preparation protocols that enable the use of quantum computers to study this question. We present protocols both for individual scar states in a particular model, as well as superpositions of them that give rise to coherent dynamics. For superpositions of scar states, we present both a system-size-linear depth unitary and a finite-depth nonunitary state preparation protocol, the latter of which uses measurement and postselection to reduce the circuit depth. For individual scarred eigenstates, we formulate an exact state preparation approach based on matrix product states that yields quasipolynomial-depth circuits, as well as a variational approach with a polynomial-depth ansatz circuit. We also provide proof of principle state-preparation demonstrations on superconducting quantum hardware.

I. INTRODUCTION
Recent decades have seen remarkable advances in our understanding of how quantum statistical mechanics can emerge from isolated strongly-interacting quantum mechanical systems. One of the most fundamental of these advances is the so-called eigenstate thermalization hypothesis (ETH) [1][2][3][4], which posits that individual quantum mechanical eigenstates at finite energy density become locally equivalent in the thermodynamic limit to equilibrium Gibbs ensembles at a corresponding temperature. Such eigenstates also govern the approach to this local equilibrium under unitary dynamics [5]. In parallel with these developments, there have been enormous strides towards realizing quantum technologies based on coherent quantum systems that are approximately isolated on experimentally relevant time scales. This gives rise to the possibility of testing the ETH and the related phenomenon of quantum information scrambling experimentally using analog quantum simulators [6][7][8][9][10] and digital quantum computers [11]. * iadecola@iastate.edu Along with this progress in understanding the ETH and quantum thermalization has come the realization that there are quantum systems that do not thermalize under certain conditions. Two notable examples include integrable [12,13] and many-body localized systems [14][15][16], where an extensive number of conserved quantities preclude the possibility of reaching a conventional locally thermalized state. An alternative means of avoiding thermalization is provided by quantum many-body scars (QMBS), a phenomenon whereby nonintegrable quantum systems exhibit a set of rare finite-energy-density eigenstates that do not satisfy the ETH [17][18][19]. Such eigenstates have been found in a variety of contexts, including the Affleck-Kennedy-Lieb-Tasaki (AKLT) model [20,21], ensembles of Rydberg atoms [22][23][24][25], and various other interacting spin [26][27][28][29][30][31][32][33][34], bosonic [35], and fermionic [36][37][38] models. These eigenstates can also give rise to coherent periodic dynamics from certain initial states, which has allowed QMBS to be observed in quantum simulation experiments [22,25,39,40].
Substantial effort has been devoted to formulating mathematical criteria for the emergence of QMBS [28,29,31,34,35,41], and all examples so far require some degree of fine tuning. Understanding the fate of scarred eigenstates and the associated coherent dynamics un-der perturbations is thus an important research direction, but relatively little progress has been made so far. Ref. [42] found a general lower bound on the thermalization timescale for a scarred state in the presence of a perturbation, t * = O( −1/(d+1) ), where is the perturbation strength and d is the system's spatial dimension. However, this bound relies only on the underlying Hamiltonian's spatial locality, and numerical studies of specific models in, e.g., Refs. [34,42] found lifetimes that substantially exceed this bound.
Studying the lifetime of QMBS under perturbations is challenging because, although scarred eigenstates often have modest entanglement and can be represented efficiently in terms of matrix product states (MPSs) [21,30,43], perturbations couple them to states nearby in energy which typically have extensive volume-law entanglement. Analytical perturbation theory is impractical here owing both to the complexity of these highly excited eigenstates and to the exponentially large density of states at finite energy density. Numerical methods to directly evaluate the real time evolution of a scarred state under perturbations also encounter challenges. Exact methods are limited to small system sizes, while approximate tensor network methods [44,45] are generally limited to early times [46].
A natural question is whether one could exploit quantum computers to investigate the behavior of scarred states under perturbations. It has long been known that quantum computers can evaluate real-time dynamics efficiently [47]; indeed, simulating quantum dynamics is one of the leading candidates for near-term practical quantum advantage [48,49]. However, even if we assume access to a noiseless quantum computer that can perform accurate time evolution, we are still left with the challenge of state preparation: how can we prepare scarred states on a digital quantum computer, and what are the resources required to do so? This is the subject of the present work.
There are two general state preparation tasks that one faces in this context, depending on whether one wants to investigate the lifetime of scarred dynamics or scarred eigenstates. In the first case, the aim is to time-evolve a superposition of scarred states that exhibits periodic dynamics in the unperturbed limit and extract the lifetime of the observed oscillations. In many cases of interest, a product state is sufficient for these purposes and the state preparation is therefore trivial [22,23,25,26,34,35,39,40,50]. However, in other cases, the simplest superposition of scar states is area-law entangled and has a nontrivial MPS representation with a finite correlation length [27,51,52], and here some thought must be put into the most efficient method to prepare such states.
In the second case, the aim is to prepare an individual scarred eigenstate and evolve it under the perturbed Hamiltonian. In many cases of interest, scarred eigenstates have entanglement scaling logarithmically with system size [21,24,26,27,36,37,51,53], similar to critical states in one dimension that are described by conformal field theories. Although the entanglement content of these states is modest compared to typical volume-law states at the same energy density, it is an open question what are the minimal quantum resources needed to prepare such a state.
In this work we address both of the above cases for a specific model with QMBS. In Sec. II we define the model and state preparation tasks in more detail. In Sec. III we consider the problem of preparing a particular class of superpositions of scarred eigenstates that can be realized as a one-parameter family of MPSs with bond dimension χ = 2. We consider two related approaches. First, we explicitly construct a linear-depth circuit that prepares the desired state with perfect fidelity. Second, we discuss a probabilistic method that prepares the desired states in constant depth using measurements and postselection. The latter method has a postselection success probability that decays exponentially with system size, albeit with a base that can be tuned by adjusting the circuit depth. This allows for a flexible tradeoff between circuit depth and success probability that is advantageous for implementation on near-term quantum hardware. In Sec. IV, we discuss two strategies for the preparation of individual scarred eigenstates. In the first strategy, we identify MPS representations for the scarred eigenstates and convert these to quantum circuits that prepare the states with perfect fidelity in quasipolynomial depth. In the second, we propose a polynomial-depth variational ansatz that we show numerically captures the scarred eigenstates with at least 99% fidelity at numerically accessible system sizes. In both Secs. III and IV, we provide proof-of-concept demonstrations of the respective state preparation tasks on Rigetti quantum processing units (QPUs), often with some additional simplifications in order to obtain better results on hardware. Finally, we provide a conclusion and outlook in Sec. V.

II. MODEL AND STATE PREPARATION TASKS
In this work, we select a particular reference model with QMBS to exemplify our state preparation techniques for scarred eigenstates and their superpositions. We consider the spin-1/2 model defined in Ref. [27], whose Hamiltonian reads This Hamiltonian acts on a chain of N qubits with open boundary conditions, and each qubit is equipped with Pauli operators X i , Y i , and Z i . Note that H 0 commutes with the operators Z 1 and Z N -the Z-basis projections of the edge qubits are therefore conserved quantities. H 0 also conserves the number of Ising domain walls, measured by the operator n DW = The λ term in Eq. (2.1) can be viewed as a kinetic term for domain walls, while the Ising interaction J serves as a chemical potential for the domain walls. The ∆ term induces nonlocal interactions between domain walls and makes the model nonintegrable. It is interesting to note that this model is dual to a Z 2 lattice gauge theory coupled to fermionic matter [54] and can be realized in Rydberg atom quantum simulators in the antiblockade regime [55].
The model (2.1) has two towers of QMBS states related by the global spin-flip operation G = N i=1 X i . The first tower is given by . . 0 , and k = 0, . . . , L/2 − 1 (we take N even for simplicity). The raising operator for this tower of states is given by has energies E k = −∆N + J(N − 1) + (2∆ − 4J)k. Both towers of states therefore have an extensive energy bandwidth, so that typical states in each tower correspond to highly excited states of the model H 0 . Nevertheless, states in either tower for which k/N is finite as N → ∞ have bipartite entanglement entropy scaling as ln(N ), in contrast to the volume-law entanglement entropy of typical eigenstates at the same energy density [27]. In the remainder of the paper we will restrict our attention to the tower {|S k }, but all results we obtain for this tower hold equally well for the tower {|S k }. Dynamical consequences of these scar states can be observed by evolving the system from a suitable family of initial states with support on the towers of interest. In Ref. [27], it was shown that the following family of states parameterized by ξ ∈ C is exclusively supported on the tower {|S k }: where the projection operator excludes any computational basis state in Eq. (2.5) containing the local configuration |1 i |1 i+1 , and where the normalization factor The state |ξ can be decomposed onto the tower {|S k } as follows: As such, for any value of ξ, evolving the initial state |ξ under H 0 yields perfect periodic revivals of |ξ with period T = π/(∆+2J) set by the energy spacing between consecutive states in the tower. Unlike typical states in the tower {|S k }, the state |ξ is area-law entangled. In fact, it can be written as an MPS with bond dimension χ = 2 [27], so its bipartite entanglement entropy is upper bounded by ln 2.
The state |ξ has connections to several interesting problems in condensed matter and atomic physics. For instance, it is unitarily equivalent to a state found in Ref. [56] to be a good approximation to the ground state of a system of Rydberg atoms in the so-called Rydberg blockade regime [57,58], where the computational basis states |0 i and |1 i correspond to atom i being in its ground state or a highly excited Rydberg state, respectively. Due to strong interactions, such systems are subjected to an energetic penalty for having two excited atoms next to one another. This constraint, sometimes called the "Fibonacci constraint" because the number of states of m qubits that satisfy it is given by F m+2 (where F is the -th Fibonacci number), is implemented by the projector P fib in Eq. (2.5). Intriguingly, this constraint also emerges in theoretical descriptions of the ν = 1/3 Laughlin fractional quantum Hall (FQH) state in a particular quasi-1D limit [59][60][61][62]. In fact, the ground state of the system in this limit is unitarily equivalent to |ξ for a particular choice of the parameter ξ [59]. Thus, our state preparation results for |ξ will also be applicable to the seemingly disparate settings of Rydberg-atom quantum simulators and FQH liquids.
Studying the stability of scarred eigenstates and dynamics on quantum computers requires algorithms for high-fidelity preparation of the states |S k and |ξ , respectively. This paper provides a survey of approaches to both state preparation tasks and a snapshot of their feasibility with current quantum hardware. We first address the preparation of the superposition state |ξ in Sec. III, before moving onto the preparation of the scar states |S k in Sec. IV.
Before proceeding, we note that the alternating signs (−1) i appearing in Eqs. (2.3) and (2.5) can be removed by the simple unitary circuit i odd Z i . Therefore, it is useful to define the states which are now equal-amplitude and equal-sign superpositions of computational basis states. We will at times find it more convenient to work with these "tilde" states than with the original states.

III. PREPARING THE STATE |ξ
There are several possible approaches to preparing the superposition state |ξ from Eq. (2.5). In Ref. [27] it was shown that |ξ is the unique ground state of a local parent Hamiltonian with finite correlation length (see also Ref. [56]). Thus one approach is to prepare |ξ adiabatically using an appropriate parameter sweep from a zero-correlation-length paramagnetic Hamiltonian. We investigate this approach in App. A, where we find evidence that the gap of the interpolating Hamiltonian does not close with increasing N and the state |ξ is a suitable candidate for adiabatic state preparation. In practice, adiabatic state preparation on a digital quantum computer suffers from Trotter error, even for a fully gapped interpolation Hamiltonian and assuming perfect implementation of the Trotter circuit. For this reason a finite depth circuit always incurs finite error. This motivates the consideration of alternative state preparation strategies.
In Sec. III A, we demonstrate that perfect state preparation can be achieved with a circuit of depth O(N ). Sec. III B shows that a stochastic strategy [63][64][65][66][67][68] using measurements and post-selection can reduce the circuit depth to a constant at the price of an exponential postselection overhead. Finally, in Sec. III C, we benchmark both state preparation strategies on Rigetti QPUs.  In this section, we show that the state |ξ; m as well as its counterpart |ξ; m with alternating phases removed per Eqs. (2.9) can be prepared in linear depth O(m) by a unitary circuit U ξ (m) consisting of (m − 1) controlled Y -rotation gates and one Y -rotation gate. We note that Linear-depth circuit showing the preparation gate U ξ (m) to prepare |ξ; m from the zero state |Ω . This circuit consists of (m − 1) control-rotation gates and one rotation gate. The rotation angles θj are given in Eq.
a similar linear depth circuit was obtained in the FQH context in Ref. [61], but that the nonunitary approach explored in Sec. III B has not been discussed in this context.
To understand this preparation circuit, we recall that |ξ or |ξ; m is the superposition of all computational basis states excluding any local configuration |1 j−1 |1 j with the weights of each state determined by the parameter ξ. Starting with all qubits in the |0 state, the absence of |1 j−1 |1 j pairs is guaranteed under a sequence of controlled Y -rotations C 0 R Y,j−1,j (θ j ) by an angle θ j targeted on the j-th qubit that is triggered only by |0 of the (j − 1)-th qubit. For a set of properly chosen rotation angles θ j , |ξ; m can be prepared by this preparation circuit: The corresponding circuit diagram of the preparation gate U ξ (m) is shown in Fig. 1. The angles are determined as a function of ξ according to the following recursion relation: The preparation circuit Uξ(m) for the state |ξ; m is obtained by simply removing the alternating phase factors (−1) j+1 from the above definition. We will prove this recursion relation for the rotation angles and the preparation circuit in the following. An alternative definition of |ξ is given in Ref. [27] as where Z is the normalization factor in Eq. (2.7). Without loss of generality, we choose the convention where j = 2 is the rightmost term in the above product, i.e., In this convention, the (j + 1)-th qubit is always in its |0 state when applying the j-th term, and hence the projector P j+1 can be omitted in the above expression. This means that With the definition of |ξ; m = N − 2 [Eq. (3.1)], we can get rid of the two qubits in the |0 state at both ends to write where we fix P 0 ≡ 1 to be consistent with open boundary conditions. Note that in this equation, the j-th qubit of |ξ; m corresponds to the (j + 1)-th qubit of |ξ .
To facilitate the derivation, we recall the projector where P 0 = 0 for open boundary conditions. Note that to prepare the state |ξ; m , one can simply pick a different definition of x j to be x j = ξ.
The term in square brackets in the above equation is in the form of a controlled-rotation gate except that (1 + x j σ + j ) is not a unitary operator. To make it unitary, we first introduce a set of real constants φ j whose values will be determined later in this subsection. We can divide each term in the product by φ j to get We recognize that the local state |1 is only generated by the raising operator σ + j . Hence, the factor 1/φ j associated with P j−1 can be moved to be associated with σ + j−1 in the previous product term. We can then write This allows us to write (3.8) If we pick the operator associated with P j−1 is unitary such that where R Y,j is the Y -rotation gate on the j-th qubit and θ j = 2 arg(1 + ix j /φ j+1 ) is the rotation angle. Rewriting Eq. (3.8) in terms of Y -rotation gates, we have Since all terms are unitary, 1 Z is just a global phase factor not relevant for the state preparation and hence can be dropped. Note also that P 0 = 0 and P 0 = 1 and hence the j = 1 term is just a Y -rotation gate. Further rewriting the above equation using the controlled-Y rotation gates, we arrive at Eq. (3.2) with θ j given by Eq. (3.3).

B. Probabilistic Constant-Depth Circuit with Postselection
While the linear-depth circuit explored in the previous section is sufficient for state preparation, it may be challenging to implement on near-term devices for large N owing to the linear circuit depth. Here, we show that it is possible to prepare the equal-amplitude superposition state |ξ; N −2 stochastically in constant depth using measurements and postselection. The idea is to prepare k m-site blocks in the state |ξ; m , which can be achieved in depth O(m) using the circuit Uξ(m) [Eq. (3.2) with alternating phases removed in the recursive formula Eq. (3.3)]. The resulting state |ξ; m ⊗k obeys the Fibonacci constraint enforced by the projector P fib [Eq. (2.6)] within each m-site block, but adjacent blocks need not obey the constraint. To "stitch" the adjacent blocks together into a state that globally satisfies the constraint, adjacent msite blocks are each coupled to an ancilla qubit using an appropriate unitary operation. Measuring the ancilla qubit and postselecting onto an appropriate measurement outcome prepares the desired state |ξ; km , which can then be trivially converted to |ξ; km using Eq. (2.9).
To illustrate, let us set ξ = 1 for simplicity (this is not strictly necessary but does simplify the analytical expressions for the states and success probabilities). As an example, consider stitching together two copies of the two-qubit state |1; 2 . The two-qubit state is given by 11) and the state we aim to prepare is (3.12) We first initialize the system in |Ω = |0000 and apply the circuit U 1 (2) to the two consecutive two-qubit blocks to prepare the state There is exactly one configuration in the above superposition that would be projected away by P fib , namely |0110 . Therefore if we apply a Toffoli (CCNOT) gate controlled on qubits 2 and 3 and targeted on an ancilla qubit initialized in the state |0 a to the above state, we obtain (3.14) Measuring the ancilla qubit in the computational basis, we obtain the desired state |1; 4 whenever the measurement outcome is 0, which occurs with probability 8/9. The same procedure can be generalized to stitch together k copies of |1; 2 into the state |1; 2k using k − 1 ancilla qubits initialized in the |0 state. Each ancilla qubit is coupled to the array of 2k primary qubits using a Toffoli gate controlled by the neighboring qubits from consecutive two-site blocks. In this way, the state of each ancilla qubit records whether the Fibonacci constraint is violated for the pair of primary qubits to which it is coupled. The probability that the Fibonacci constraint is satisfied (such that the ancilla register is in the state |0 . . . 0 a ) is given in terms of Fibonacci numbers F by F 2k+2 /F k 4 , which decays exponentially with k. The above success probability can be improved by using the same strategy to stitch together larger blocks. For example, suppose we wish to prepare the state |1; 2m by stitching together two m-site blocks prepared in the state |1; m using the O(m)-depth circuit U1(m). The stitching can again be achieved by coupling the two neighboring qubits from the two blocks to an ancilla qubit using a Toffoli gate. An illustration of this procedure for general ξ is shown in Fig. 3. The success probability can be obtained by noting that there are F 2m+2 states of the full 2m-qubit system that satisfy the Fibonacci constraint (for which a measurement of the ancilla qubit would yield 0), while the initial tensor product state contains F 2 m+2 configurations; the success probability is then F 2m+2 /F 2 m+2 . Applying the same logic to k blocks of m sites yields a postselection success probability   k for fixed m, it grows with m at fixed k. Thus, the exponential sampling overhead can be mitigated by increasing m at the price of increasing the depth of the state preparation circuit. Fig. 2 shows the success probability for a lattice of length N using blocks of length m. On present-day NISQ devices, gate error and qubit decoherence rates are sufficiently high that the reduction of circuit depth at the expense of postselection may be an acceptable tradeoff. We analyze this tradeoff further in Sec. III C when we implement this state preparation protocol on quantum hardware.

C. QPU Results
We now implement the |ξ state preparation protocols on Rigetti's Aspen QPUs and discuss the execution results of the linear-depth circuit in Sec. III A and its probabilistic post-selection variant in Sec. III B. We set ξ = 1 for numerical evaluations, in which case the Pauli rotation angles in Fig. 1 can be summarized as The quilc compiler generates an optimized program for the instruction set architecture of QPU chips, compiling all logical gates into the following group of native gates: RZ(θ), RX(π/2), RX(π), CPHASE(θ), CZ, XY(θ).
Results obtained from the NISQ hardware are typically prone to multiple sources of error, leading to deviations from ideal unitary calculations. Here we quantify these errors by estimating two figures of merit: first, the Bhattacharyya distance that captures the divergence between the exact and the observed distributions, p(x) and q(x), of the measured bitstrings x. Second, we use the expectation value of the parent Hamiltonian of the state |ξ , namely [27] that vanishes H ξ = 0 in the case of ideal noiseless preparation of the state |ξ [69]. We collect 10 observed samples of these quantities, along with the success probability of the post-selection protocol, for a range of different system sizes N and block sizes m. Each sample measurement is made with 10 4 shots for the energy and the success probability, and 10 5 shots for the Bhattacharyya distance. The individual samples are depicted as round dots in Fig. 5 and 6, and their average values are connected with solid lines. See the captions of Fig. 5 and 6 for the specification of used device nodes. When the qubits are measured immediately after state preparation, without performing further unitary operations for, e.g., time evolution, the two-step process of applying Toffoli gates and post-selecting on ancilla bits can be reduced to the classical post-selection of N -qubit output bitstrings. Such replacement reduces the room for error, since Toffoli (CCNOT) is a non-native gate that must be compiled into multiple imperfect two-qubit gates (e.g., CPHASE, CZ, XY) before running on QPUs [70]. It also "unstitches" the whole circuit back into N −2 m decoupled blocks of U ξ (m), which introduces certain computational advantages. For example, the asynchronous execution of fragmented circuits allows us to simulate systems larger than the hardware size, or to reduce the overall error by avoiding the usage of low-fidelity qubits, at an exponential overhead of classical post-processing [71,72]. With this in mind, we repeatedly use the same m-qubit sublattice that exhibit the best average readout and gate fidelities, and then classically combine the measurement outcomes with the Fibonacci constraint to simulate N = αm + 2 qubit observables (α ∈ N). Fig. 5 displays the results of the N = 14 state preparation experiment with different block sizes m, to which the circuit U ξ (m) is applied. We first discuss the raw values without error mitigation, shown as green dots. As m increases, both the energy expectation value H ξ and the Bhattacharyya distance deviate farther from their ideal values of 0. At the same time, the success probability grows with increasing m. These manifest the trade-off between the error and the post-selection success probability: the use of a greater number of smaller circuit fragments facilitates higher-fidelity execution of the circuit, while demanding an increased sample complexity. Note that our experimental protocol is a special case of quantum divide-and-conquer [71,72], where each block does not depend on another block's measurement result. Fig. 6 summarizes the output of the 4 ≤ N ≤ 14 state preparation circuit that merges N −2 2 copies of U ξ (m = 2). As the system size N increases, the Bhattacharyya distance stays close to 0. While the energy deviation shows an apparent O(N ) scaling, its values are considerably less than those with larger m building blocks. We also illustrate in Fig. 7 the state tomography of |ξ = 1 ξ = 1| in the computational basis. In particular, its right panel highlights the sparsity of |ξ states at N = 14.
Finally, we attempt to improve the precision of experimental results by adopting error mitigation techniques, which produce the blue and orange dots in Fig. 5 and 6. To handle the measurement error, we apply the Bayesian unfolding [73] that uses the confusion matrix of 8-bit strings associated with Aspen-M's octagonal layout. The samples obtained from "corrected" bitstrings are colored in blue.
They amplify the error measured in the Bhattacharyya distance and underestimate H ξ when compared to the raw samples.
We note that, for the U ξ (m = 2) block calculation, most blue samples have negative energy that violates the positive semi-definiteness of the Hamiltonian Eq. (3.18). The presence of the spurious samples is apparently an artifact of the Bayesian unfolding. It likely occurs due to a phase error since the sample Bhattacharyya distances are close to 0. When using larger circuit blocks U ξ (m > 2), however, the Bayesian unfolding, combined with randomized compilation [74] and readout symmetrization [75], helps to lower the estimated energy H ξ while maintaining the Bhattacharyya distance near the same level. An individual orange dot represents the average measurement results from 30 logically equivalent circuits varied by random twirling gates between each cycle and ran-dom Pauli operators inserted at the circuit end along with classical post-processing that reverts the induced basis changes.

IV. PREPARING THE STATES |S k
Since generic states |S k in this tower have entanglement scaling logarithmically with system size (see Sec. II), they cannot be prepared with quantum circuits of constant depth. In this section, we demonstrate two (quasi-)polynomial-depth algorithms for preparing these states. The first, described in Sec. IV A, relies on building an MPS representation of the states |S k and then converting this representation to a quasipolynomial depth quantum circuit. The second, described in Sec. IV B, is a variational strategy that uses a polynomial-depth variational ansatz circuit to represent the states. Finally, we provide a proof-of-principle realization of these states on quantum hardware in Sec. IV C, where we also write down a simplified linear-depth circuit for preparing the highest-weight state in the tower, |S N/2−1 .

A. Quasipolynomial Depth Circuit from Matrix Product State Representation
To facilitate the discussion below, we rewrite the tower of scarred eigenstates for N sites as follows [see Eqs.  where | s = m i=1 |s i are computational basis states labeled by bitstrings, M si are square matrices of size χ, and L| and |R are respectively χ-dimensional row and column vectors implementing open boundary conditions. One can view an MPS as a representation of a deterministic finite automaton (DFA), where the matrices M si correspond to a transition matrix M for the DFA states [79]. Since, for any i, M si is a square matrix of size χ, this DFA will have χ states, and M si j,k will be nonzero if the character s i takes the state j to the state k. The bra L| and ket |R denote the initial and final states of the DFA, respectively. It is important to note that the states in the DFA live in the auxiliary bond space of the MPS and are not quantum states themselves.
DFAs are used to represent regular expressions, or string expressions using a finite set of characters. To see how we can use this language to represent quantum states, consider |D 4 2 : The set of computational basis states in the superposition forms a language {|0101 , |1001 , |1010 } from which a regular expression can be constructed, where the characters in the regular expression are spin states {|0 , |1 }. For the states |D * k = ∞ n=0 |D n k , the corresponding regular expression is where * is a Kleene star (a * = I + a + aa + aaa + ... where I is identity (or empty string)) and is only applied to |0 , and [.] m means to repeat the term inside the bracket m times.
An example of a DFA representing the regular expression for |D * k (with k = 4) is shown in Fig. 8. Note that if we limit the number of transitions through the DFA to m, the regular expression for |D * k reduces to that for |D m k . For general k, the DFA is defined on 2k states {S 0 , A j , B j , F k } j=1...k−1 . The DFA starts at S 0 and ends at F k . The first computational state can either be |0 , which transitions the DFA back to S 0 , or |1 , which transitions the DFA to A 1 . We denote this by For the intermediate states, we have Finally, the DFA terminates at the final state F k , where All these transition rules can be summarized by a transition matrix M, where M i,j gives the character needed to take the DFA from state i to state j: The full transition matrix is then To obtain the matrices from Eq.
That is, M 0 is obtained by applying 0| to every entry in M, and M 1 is obtained by applying 1| to every entry in M. These are 2k × 2k matrices, and so the overall bond dimension of the MPS will be 2k. Also, since we start at state S 0 and end at state F k , the boundary vectors become L| = 1 0 · · · 0 |R = 0 · · · 0 1 T (4.11) One can now read Eq. (4.2) from left to right as a traversal through the DFA in Fig. 8, starting at DFA state S 0 , making m transitions according to the spin configuration s, and ending at DFA state F k .
The resulting MPS can be converted into a quantum circuit by following the sequential preparation procedure of Ref. [77]. Orthogonalizing the MPS in one direction results in m two-level unitary matrices of width log 2 4k . Using Gray codes [80], each of these unitaries can be decomposed into O(k) controlled-unitary gates, each spanning over O(ln k) qubits [81]. Each of these controlledunitaries then requires O(ln 2 k) CNOT gates. The depth of the resulting circuit is then O mk ln 2 k . Fig. 9 depicts an example preparation circuit for the state |D 6 2 . The circuit makes use of the following single-qubit uni-taries: (4.12)

B. Polynomial Depth Variational Ansatz
A strategy similar to the recursive linear circuit used in the preparation of the |ξ state suggests a variational circuit architecture for creating the |S k state should exist. Up to unimportant relative phase factors, the |S k state is an equal-weight superposition of all bitstrings of N sites, such that there are k 1s, the first and last sites are 0, and there are no two neighboring 1s. All such  basis states can be generated by a four-qubit buildingblock unitary operator that transforms the state |0100 to a linear superposition of states |0100 and |0010 [82], and acts trivially on all other 14 basis states. One can construct the unitary by applying a two-qubit gate on the two middle sites controlled by the first and the last qubit such that the gates act nontrivially only when the first and the last qubits are 0, in which case the two middle qubits transform as We have chosen a building block with only real elements for easier optimization. Mathematically, we can write the unitary U j (θ j ) acting on qubits j to j + 3 as as (4.13) where the operators P j = |0 j 0| j implement the control on the first and last qubits and (X j+1 Y j+2 − Y j+1 X j+2 ) acts as a Pauli-Y on the 01 and 10 states of the two middle qubits, while annihilating 00 and 11. The building block U j is illustrated in Fig. 10(a).
We now construct an ansatz circuit. Since the gates are designed to conserve both the total number of 1s, i.e., k, and the Fibonacci constraint of no neighboring 1s, we start with the following initial state that has k 1s: (4.14) In Fig. 10(b), we show an example of the circuit acting on the above initial state for k = 2 and N = 8. Our ansatz consists of multiple staircases of the U j gates. In the first layer, we apply U 1 ...U 2k−3 U 2k−1 to move the last 1 from qubit 2k to qubit 2k + 1 and generate all states with the last 1 at or before qubit 2k + 1. The next layer, U 2 ...U 2k−1 U 2k , generates all states with the last 1 at or before 2k + 2. These staircases continue until we reach the staircase starting with U N −3 , acting on the last four qubits. All these staircases have unitaries acting on every other site. These layers generate all d states. However, we have found that the number of variational parameters is insufficient to obtain an equal-weight superposition. Thus the ansatz also contains a complete staircase at the end with a unitary U 1 U 2 ...U N −3 to equalize the probabilities. We then optimize the variational angles in these four-qubit unitaries to generate an equal-weight superposition of all d basis states by minimizing 1 − | ψ|S k | 2 , where |ψ is the state after the application of the variational unitaries. A final layer of Z gates corrects the signs of the amplitudes. Two important characteristics of this circuit are the dimension d of the constrained Hilbert space and the number n a of variational parameters. While d scales exponentially in the asymptotic limit, n a is quadratic in system size. The dimension of the constrained Hilbert space is . We need our circuit to generate all the states subject to the Fibonacci constraint and the conserved number of 1s. Thus, the basis states can be viewed as k blocks of 01 as partitions between N − 2k − 1 zeros with one zero tagged at the end of the chain. The number of variational angles n a is easy to calculate. For even N , the first and the second layer have k gates each, the third and fourth layer k + 1 gates, and so forth until we reach the staircase with N/2 − 2 gates containing every site from 2 to N − 4. We then have a staircase with N/2 − 1 gates and then a full staircase with N − 3 gates. Thus n a = 2[k+(k+1)+...+N/2−2]+(N/2−1)+(N −3), which simplifies to Similarly for odd N , we have n a = 2[k + (k + 1) + ... + (N + 1)/2 − 2] + (N − 3), which leads to Importantly, the number of variational gates is quadratic in both N and k.
The numerical optimization results are listed for several N and k in Table I. The system exhibits many local minima (local minimization routines with different initial angles tend to converge to different optimal angles), making global optimization challenging. However, we have found that constrained local optimization with random initial angles between 0 and 2π finds solutions with reasonably small error 1 − | ψ|S k | 2 < 0.01. For many cases with smaller n a and d, and less complex optimization, the minimum error is zero to within machine precision, suggesting the |S k state may be exactly reachable with the ansatz architecture. In practice, we can only find very good approximations for large d and n a .

C. QPU Results
It is noticeable that the doubly-controlled gates are essential pieces in building up the state preparation circuits described in Secs. III A and IV B. Compiling them into the set of native gates often produces quantum programs relatively expensive in the two-qubit gate complexity. For experimental realization of the |S k state preparation, we therefore consider the specific case of k = k max ≡ N/2−1 for which the alternative circuit in Fig. 11 that involves only (N − 3) two-qubit gates can be used.
We can deterministically realize the maximal k condition by setting the spin state of the (i + 1)-th qubit to be opposite to that of the i-th qubit for all i ∈ {2, 4, · · · , N − 2}. This can be implemented by initializing the target qubit in |0 and then acting with the controlled bit-flip operator, triggered only when the control qubit is in the |0 state. Next, the superposition state can be prepared with a variant of the unitary block Eq. (3.2) where we modify the controlled rotations to be triggered only if the control qubit is in the |1 state. Running it on the even-indexed qubits between 2 ≤ i ≤ N − 2, we obtain a superposition of the following bitstring states, where c a,b ∈ R and a,b c 2 a,b = 1. Lastly, Eq. (4.19) becomes the equal-amplitude superposition state |0 ⊗ |D N −2 kmax ⊗ |0 if the Pauli rotation parameters are (4.20) Figure 11 illustrates the |S kmax state preparation circuit. We execute the state preparation circuit on Aspen M-2 and M-3. Ideally, if we perform Pauli-Z measurements after running the circuit, the resulting N -bitstrings must include k max non-consecutive 1's and begin and end with 0 by construction. However, due to limitations of current NISQ hardware, the circuit simulation of the |S kmax state often induces an error in the total spin measurement k = k max and may violate the boundary condition or the Fibonacci constraint. Therefore, to increase the fidelity of the simulation, we post-select the measurement bitstrings along with optional application of the readout error mitigation and randomized compiling techniques.
Assuring the total spin selection rule allows us to introduce the compact encoding scheme that uses only k max qubits, i.e., indexed as {2, 4, · · · , N −2}, to encode the Nqubit states satisfying k = k max without consecutive 1's. Before After Before After Before After Before After It is equivalent to decoupling the non-entangled boundary spins and replacing the error-prone quantum operator Eq. (4.17) with deterministic post-processing. Its use of fewer qubits also brings the advantage of circumventing non-linear qubit connectivity and considerable error reduction. To make sure the Fibonacci constraint is satisfied, we apply the following projection operator that removes |01 configurations in the compressed k maxqubit encoding, or equivalently, |0110 configurations in the full N -qubit representation.
Having run the state preparation program, we measure the errors through the Bhattacharyya distance defined in Eq. (3.17) between the ideal and empirical bitstring distributions and the expectation value of the spin-1/2 Hamiltonian from Eq. (2.1). With the k max -qubit encoding, the Hamiltonian is mapped to where H λ = 0, H ∆ = 2, and It is obtained by applying the encoding rules of Table II, where vanishing substitutions on the second and third columns impose the total spin selection rule k = k max . This reduces the task of computing H 0 to measuring H J . Fig. 12 collects the results of the |S kmax state preparation experiment, showing the Bhattacharyya distance, the energy H J and its error, and the success probability of the Fibonacci projection, see Eq. (4.21). Individual sample points are drawn as round dots, and their averages are connected as the solid line. Each dot aggregates 10 4 shots for the energy and success probability, and 10 5 shots for the Bhattacharyya distance. We record the qubit indices used for the QPU experiment in the caption of Fig. 12.
We consider the raw samples without the readout error mitigation and randomized compiling techniques, denoted as green dots. We note a rapid drop of the success probability as N increases. Unlike the probabilistic |ξ state preparation protocol in Section III B, the noiseless probability should be 1, and thus the decrease from 1 highlights the detrimental effects of the noise. While the measured Bhattacharyya distance after the Fibonacci projection remains close to 0, the estimated energy H J exhibits a linearly growing error rate as compared to the ideal value −(N − 3), reaching around 10% at N = 14. Also note that the noise can create spurious samples whose H J is less than the theoretical minimum of the H J in Eq. (4.23) Next, we try various error mitigation techniques to improve the estimate accuracy. To reduce the measurement error, the iterative Bayesian unfolding [73] uses the confusion matrix of 8-bit strings that reflects Aspen-M's octagonal layout. The samples obtained from the modified bitstrings are marked as blue dots and feature interesting changes. They have a higher success probability-i.e., they remain more in the Fibonacci subspace-and feature less error in estimating H J after the Fibonacci projection. However, their Bhattacharyya distance also amplifies significantly with increased sample variance. Application of randomized compilation [74] and readout symmetrization [75] can help to reduce the Bhattacharyya distance while maintaining other improvements, as seen from the orange entries in Fig. 12. Figure 13 displays the results of the |S kmax state tomography after applying the aforementioned error mitigation methods and the Fibonacci post-projection. It illuminates the sparsity of the |S kmax state that populates only N/2 bitstrings out of 2 kmax different choices.

V. CONCLUSION AND OUTLOOK
In this work we have explored various approaches to preparing quantum many-body scar states and their superpositions on quantum computers with a focus on the spin-1/2 chain model of Ref. [27], where an emergent Fibonacci constraint ensures that superpositions of scarred eigenstates must be entangled. Linear-depth circuits to prepare a one-parameter family of area-lawentangled superpositions of scarred eigenstates, as well as a nonunitary preparation scheme that uses measurement and postselection to prepare the same superposition in constant depth are provided. We also derived an MPS representation of the tower of scarred eigenstates and used this to generate a quasipolynomial-depth circuit to prepare individual scar states; additionally, we proposed a variational scheme based on a polynomial-depth ansatz that captures the scar states with fidelity at least 99% in all numerically accessible cases. Proofof-concept demonstrations of scar-state preparation were executed on quantum hardware.
While the state preparation protocols explored in this work were formulated with a specific model in mind, they can readily be adapted to other towers of scar states. For instance, entangled superposition states analogous to |ξ are known for both the AKLT model [52] and the bondbimagnon tower in the spin-1 XY model [26,51]. In both cases, the initial state is a finite-bond-dimension MPS which can be viewed as the projection of a simple wavefunction on χ-dimensional qudits. One expects that nonunitary methods along the lines of Sec. III could be used to prepare such states. Likewise, both towers admit exact MPS representations that can be converted to quasipolynomial-depth quantum circuits. It is also worth noting that simpler state preparation protocols can be used in some cases-for example, the bimagnon tower in the spin-1 XY model [26] can be superposed into a product state that is trivial to prepare, and the scarred eigenstates themselves are effective spin-1/2 Dicke states that can be generated using a slight variation on the approach of Ref. [76].
Several state preparation techniques are worth future investigation. First, for the superposition state |ξ , it is desirable to find a deterministic (i.e., postselectionfree) finite-depth state preparation protocol in which unwanted ancilla measurement outcomes are corrected by local unitary operations instead of being discarded. Such a scheme can leverage the fact that the local form of the state in the vicinity of an ancilla measurement outcome of 1 is fixed. For example, consider stitching together two states of the form |1; m with a single ancilla measurement, and suppose the measurement outcome was 1. Then the primary qubit register is in the state since the state away from the central two sites (which are projected onto |11 by the measurement) obeys the Fibonacci constraint. One can imagine trying to correct the four site block |0110 using a unitary circuit controlled by the states of the two adjacent qubits from the (m−2)-site blocks to its left and right. However, the target states in the four cases (labeled by the four possible states of the adjacent qubits) have different normalization factors due to the Fibonacci constraint; therefore, the state of the qubit chain after feedback will not be an equal-amplitude superposition state as desired. For example, this situation should be contrasted with a recently proposed finite-depth deterministic scheme to prepare the AKLT ground state [68], which is able to correct undesired measurement outcomes with a local circuit. However, the correction operation designed in that work makes use of the fact that the AKLT state is a symmetry-protected topological state [83,84], which is a property that is not shared by the state |ξ . We leave the question of whether an undesired measurement outcome in our probabilistic protocol can be corrected by a finite-depth circuit for future work. With regard to the states |S k , it would be interesting to determine whether the variational ansatz circuit explored in Sec. IV B can represent the |S k states exactly. This would conclusively demonstrate that the scarred eigenstates of the model 2.1 can be prepared in polynomial depth (as opposed to quasipolynomial depth as shown in Sec. IV A), which is known to be possible for, e.g., Dicke states [76]. Attempts to find an exact solu-tion for the rotation angles in the circuit architecture of Sec. IV B yield n a coupled nonlinear equations that are intractable in all but the simplest cases. It would be interesting to see whether a modified circuit architecture involving gates of the form in Eq. (4.13) can be used to obtain an analytically tractable system of equations for the rotation angles.
It may also be worth considering whether the recursive method of Ref. [76] to prepare Dicke states can be adapted to exactly prepare the |S k states, or equivalently the projected Dicke states |D m k [see Eq. (4.1)], in polynomial depth. The recursive construction of a state preparation circuit in Ref. [76] hinges on a recursion relation for Dicke states with different magnetizations. The projected Dicke states also obey a recursion relation: However, unlike the recursion relation for Dicke states, this one relates |D m k to both |D m−1 k and |D m−2 k−1 , which are defined on qubit registers of different sizes. This complicates the possibility of a straightforward modification of the protocol of Ref. [76], and we leave this question for future work.
Another potentially useful avenue to explore is the application of nonunitary methods to prepare the scarred eigenstates |S k . In particular, these eigenstates can be prepared from the state |ξ by projectively measuring the magnetization operator M z = N i=1 Z i . For any k = 0, . . . , N/2−1, |S k is an eigenstate of M z with eigenvalue N − 2k. Since each of these eigenvalues is unique, we see from Eq. (2.8) that the projection of |ξ into an M z eigenspace must yield one of the |S k . The probability of a given measurement outcome k is peaked around a value that depends on |ξ| 2 ; this most probable outcome can be tuned from k = 0 as |ξ| 2 → 0 to k = N/2 − 1 as |ξ| 2 → ∞. Thus, given the ability to measure M z without fully collapsing the system into a z-basis eigenstate, one can in principle prepare the state |S k in constant depth. One way to achieve such a measurement is by a dispersive coupling χ d ( i Z i )a † a between the qubits and a single cavity mode, where the measurement outcome would be recorded as a shift of the cavity frequency by an integer multiple of χ d . A projective measurement of M z can also be achieved using an ancilla register containing log 2 N qubits, one for each digit of the binary representation of the magnetization eigenvalue, using methods proposed in Refs. [85] and [86]. This state preparation method also suffers from postselection overhead beyond what is needed to prepare |ξ : the probability of the most probable measurement outcome decays as a power law in N for large N . However it may have the advantage of comparatively modest circuit depth relative to the unitary state preparation protocols considered here.
Our results can feed into future studies of the stability of scarred eigenstates and their dynamics on quantum computers. There are two experiments one would like to perform-one in which the state |ξ is evolved under a perturbation of the Hamiltonian in Eq. (2.1), and another in which the state |S k is evolved. In the former case, the goal is to extract the lifetime of the oscillatory scarred dynamics from the time series of a local operator's expectation value [87], while in the latter it is to extract the lifetime of the eigenstate from similar data. One relevant class of perturbations to consider is an external magnetic field in the x-direction, i X i , which breaks the conservation of the Ising domain wall number n DW by disrupting the balance of X and ZXZ in the first line of Eq. (2.1) [27]. The evolution operator over a small time step δt can be written in a Trotter product form as (setting J = 0 for simplicity) U (δt) ≈ e −i∆δt i Zi e iλδt i odd Zi−1XiZi+1 e −i(λ+ )δt i Xi e iλδt i even Zi−1XiZi+1 . (5. 3) The nontrivial ZXZ rotations above can be compiled using standard methods into four CNOT gates, two H gates, and an R Z gate. An informative figure of merit to track while evolving the state |ξ is the expectation value of, e.g., Z i−1 X i Z i+1 , which exhibits coherent oscillations under evolution by H 0 that should decay in the presence of a nonzero . For |S k , it is informative to track the expectation value of M z , which is not conserved under H 0 for generic initial states but is conserved for the initial state |S k ; thus tracking M z (t) provides information about the fidelity decay of |S k . These calculations are in principle straightforward to carry out on present-day quantum computers, but our QPU results in this work demonstrate that a significant fraction of the QPU's coherence budget will be expended on state preparation. Thus, any near-term demonstration of this calculation must employ substantial problem-tailored optimizations and error mitigation strategies, which should be explored in future work.  Results are obtained using exact diagonalization with a total number of steps ns = 10 3 and δs = T * /ns, where T * is the minimal T value to reach at least 99% fidelity. and the target Hamiltonian H P , Note that H P is simply H ξ=1 /2 [see Eq. (3.18)], plus terms involving the projectors P i = |1 i 1| i that ensure the ground state is in the Fibonacci Hilbert space. Both Hamiltonians are positive semidefinite (i.e., their ground states have zero energy) and conserve the Z-basis projection of the first and last qubit, which can be fixed to be in the 0 state. The ground state of H X is then the product state |ψ(0) = |0 ⊗ |+ − + − · · · + − ⊗ |0 , where |± satisfy (1 ∓ X i ) |± = 0, while the ground state of H P is |ξ = 1 . To adiabatically prepare the desired state, we evolve |ψ(0) under the time-dependent interpolating Hamiltonian such that H(0) = H X and H(T ) = H P . Provided H(s) has an energy gap for all s, then if T is sufficiently large, the time-evolved state |ψ(T ) ≈ |ξ = 1 to a good approximation. The adiabatic time evolution was carried out with exact diagonalization of the Hamiltonian (A3). After |ψ(0) is prepared, the operator U (s, s + δs) is applied n s = T /(δs) times with s = nδs, 0 ≤ n ≤ n s − 1.
Here, n s is the number of iterations, T is the total evolution time and δs is the time separation. The time evolution operator U (s, s + δs) = e −i(δs)H(s) and the overlap | ψ(s)|ξ = 1 | are calculated at each iteration. In Fig. 14, we plot the difference between the instantaneous energies of the ground and first excited states E 1 −E 0 of H(s) against the interpolation parameter, s/T , for several values of N . The minimum energy gap occurs at s = T and appears to saturate to a constant as N increases. In fact, an exact calculation of the energy gap of H P gives 1 − 1/ √ 2 = 0.292 . . . [56], which appears to be in reasonable agreement with the limiting value of the numerically observed gap as obtained from a quadratic extrapolation in 1/N (see Fig. 14 inset). This indicates that high-fidelity adiabatic state preparation should be possible in a finite total time T . In Fig. 15, we plot the total time T * required to achieve state preparation with 99% fidelity as a function of N . T * is calculated by solving for the root of the function f (T ) = | ψ(T )|ξ = 1 | − 0.99 using Newton's method. We find that T * appears to approach a constant value T * ≈ 25. This is consistent with the behavior of the gap observed in Fig. 14.
These results lead to the conclusion that adiabatic state preparation should be possible for the state |ξ = 1 . However, it is important to compare the quantum resources required to implement adiabatic evolution on a QPU against those of the methods discussed in Sec. III. In particular, if the evolution operator is approximated using a Trotter product formula, the number of Trotter steps required to achieve Trotter error after evolution by a time t scales as t 2 / . Inserting t = T * ≈ 25 and demanding a Trotter error = 0.01 yields an estimated 62 500 Trotter steps. Thus the needed circuit depth is constant, but prohibitively large for execution on presentday QPUs. This fact is a well-known shortcoming of the adiabatic algorithm [88] which partially motivated the development of alternative approaches like QAOA [89].