Monitoring Quantum Simulators via Quantum Non-Demolition Couplings to Atomic Clock Qubits

We discuss monitoring the time evolution of an analog quantum simulator via a quantum non-demolition (QND) coupling to an auxiliary `clock' qubit. The QND variable of interest is the `energy' of the quantum many-body system, represented by the Hamiltonian of the quantum simulator. We describe a physical implementation of the underlying QND Hamiltonian for Rydberg atoms trapped in tweezer arrays using laser dressing schemes for a broad class of spin models. As an application, we discuss a quantum protocol for measuring the spectral form factor of quantum many-body systems, where the aim is to identify signatures of ergodic vs. non-ergodic dynamics, which we illustrate for disordered 1D Heisenberg and Floquet spin models on Rydberg platforms. Our results also provide the physical ingredients for running quantum phase estimation protocols for measurement of energies, and preparation of energy eigenstates for a specified spectral resolution on an analog quantum simulator.


I. INTRODUCTION
In the standard scenario of analog quantum simulation, a broad and tunable class of many-body Hamiltonians of interest is designed based on the resources provided by a particular physical platform. Examples in different physical platforms include spin models realized with Rydberg tweezer arrays [1][2][3][4][5][6][7] (for review see [8]), trapped ions [9][10][11] or superconducting qubits [12,13], or Hubbard models realized with bosonic and fermionic atoms in optical lattices [14][15][16]. The physical realization of these Hamiltonians then allows the study of equilibrium and non-equilibrium phenomena, where the quantities of interest characterizing the quantum many-body state are spin-or site-occupation correlation functions. These are inferred from (destructive) site-resolved readout of spins, or in the case of Hubbard models from quantum gas microscopy.
In contrast, we will be interested below in a setting where we learn the state and dynamics of the manybody system by entangling the state of the quantum simulator with an auxiliary quantum system, followed by measurement of the auxiliary degrees of freedom, as illustrated in Fig. 1(a). In its simplest form, this auxiliary quantum system is a single qubit acting as a 'clock' [17,18] (or 'control') qubit, which can be manipulated by single qubit operations (rotations). However, the following considerations generalize immediately also to a multiqubit setting, where the auxiliary system represents a small scale quantum memory or quantum computer.
At the heart of our considerations is the quantum nondemolition (QND) Hamiltonian (1) * Equal contribution which generates the QND gate entangling the quantum simulator with a clock qubit. To be specific, we assume a spin model with Hamiltonian H spin for the simulator, and we denote by {|0 c , |1 c } the logical states of the clock qubit [19]. The above Hamiltonian is QND, with H spin the 'energy' of the quantum many-body system, which plays the role of the QND variable. To illustrate the action of the above QND gate U QND (t), consider a quantum simulator prepared in superposition |ψ spin = c | of energy eigenstates, H spin | = E | with eigenenergies E . Under the QND gate, an initial state |Ψ(t = 0) = |ψ spin ⊗ |+ c of the joint quantum simulator prepared in |ψ spin and the control qubit prepared in the superposition state |+ c = 1 √ 2 (|0 c + |1 c ) will evolve according to Thus the superposition of many-body energy eigenstates gets entangled with the phase of the Bloch vector of the clock qubit rotating in the xy-plane on the Bloch sphere. A readout of this phase via the clock qubit provides us with a QND measurement of 'the energy' of the quantum many-body system. In a broader context, we note that this QND gate is also the basic building block of quantum algorithms like quantum phase estimation (QPE) [20][21][22] to measure energies of the many-body spin system, and prepare corresponding eigenstates with a prescribed spectral resolution. We will use such features below in a protocol to measure the spectral form factor (SFF), where we access correlations between eigenenergies E , encoded through a Fourier transform, via the clock qubit.
In the present paper we will first address the challenge of implementing the QND Hamiltonian (1)  for a broad class of freely designable spin models in a quantum simulator in atomic physics setups. Remarkably, as we show in Sec. II, this can be achieved with Rydberg tweezer platforms [1][2][3][4][5][6][7] using laser dressing schemes [23][24][25][26][27][28], and by employing a Rydberg blockade mechanism between the clock qubit and the simulator atoms to implement the QND Hamiltonian (1). Second, we wish to explore and illustrate the application of quantum protocols, which build on the above QND gate, providing access to novel quantum many-body observables of interest under experimentally realistic conditions. A relevant example is provided by measurement of the spectral form factor (SFF), as discussed in context of [29][30][31][32], where the aim is to identify signatures of ergodic vs. non-ergodic dynamics without an explicit spectroscopic study [33][34][35] of the energy spectrum. In Sec. III we describe a protocol, where the SFF can be measured via QND couplings to a clock qubit. We illustrate this protocol and its performance with simulated measurement runs for the disordered 1D Heisenberg model and Floquet models, which can be implemented with our techniques on Rydberg platforms.
The paper is organized as follows. In Sec. II we discuss implementation of the QND Hamiltonian with Rydberg tweezer arrays. The properties of the SFF in Hamiltonian systems are briefly summarized in Sec. III, where we also present the protocol for preparation of the initial state (Sec. III B 1) and the protocol for the SFF measurement together with the discussion of its experimental limitations (Sec. III B 2). The SFF of Floquet systems and the corresponding measurement protocol are considered in Sec. IV, and we conclude in Sec. V.

II. PHYSICAL REALIZATION OF HQND IN RYDBERG TWEEZER ARRAYS
The challenge is to implement H QND , Eq. (1), for a broad and tunable class of spin models H spin on a given physical platform. The relevant example below is the disordered 1D Heisenberg spin-1/2 model, or, to be more specific with designable single particle and two-body (interaction) terms. In implementing H QND we are required to implement the two and three-body terms involving the control qubit as H QND = H spin ⊗ |0 c 0|. Rydberg tweezer arrays in 1D, 2D and 3D [1][2][3][4][5][6][7] provide the tools to design a broad class of spin Hamiltonians H spin via Rydberg dressing [23][24][25][26][27][28]. Here long-lived atomic (hyperfine) ground states are trapped in the optical tweezer, and play the role of the spins in our quantum simulator, while the two-body interactions are engineered by admixing weakly to the ground state via off-resonant laser dressing the (strong) van der Waals interactions between Rydberg states [36,37]. Remarkably, the same 'dressing toolbox' which allow to design H spin also provide us with a recipe to engineer H QND .
As originally discussed in [38], see also [39][40][41][42], interactions engineered via laser dressing can be turned on and off by preparing a control spin (or qubit) which is in the ground state |0 c or a Rydberg state |1 c , respectively [see Fig. 1(b) left and right panel, respectively]. If the control spin is in the ground state, the control spin does not interact with the system spins, and thus the Hamiltonian H spin is realized, as discussed above. On the other hand, for the control spin in the Rydberg state |1 c , the long-range character and strength of Rydberg-Rydberg interactions will, via the dipole blockade mechanism, detune the Rydberg states of the simulator atoms, thus effectively turning off the dressing interactions, i.e. we have H QND = H spin ⊗ |0 c 0|. We will analyze this below in a realistic atomic physics setting.

A. Hamiltonian for the simulator and control qubit
We consider the setup outlined in Fig. 1(b). L atoms trapped in optical tweezers are arranged in a ring repre-senting a 1D quantum simulator of spin-1/2 with periodic boundary conditions. The distance r between simulator atoms is assumed to be larger than 2.4 µm as discussed in Appendix (A 3). The control qubit is represented by an atom trapped in the center of the ring. In Fig. 1(b) we show the atomic level structure for the atoms representing the quantum simulator, and the control atom. The spin-1/2 degrees of freedom of the simulator are encoded in two long-lived hyperfine ground states |g ± with energies E g± . These are coupled to Rydberg states |g ± → |r ± with energies E r± by two off-resonant lasers with respective frequencies ω ± and detunings ∆ ± = ω ± − E r± + E g± Ω ± , and Ω ± corresponding Rabi frequencies. The ground state of the control qubit is |0 c with energy E gc , while |1 c is a Rydberg state with energy energy E rc . We drive the control atom transition with a laser of frequency ω c which is tuned near resonance with detuning ∆ c = ω c − E cr + E gc . The corresponding Rabi frequency is Ω c . In our protocol, the control qubit is prepared initially in a superposition state of ground and excited state with a π/2-pulse.
The Hamiltonian of the total system is H tot = H s + H c + H sc , written as sum of the simulator and control atom Hamiltonians, and the simulator-control interaction. These are given by Here H vdW is the van der Waals interaction between the Rydberg manifolds {|r ± i } and {|r ± j } of simulator atoms i and j, and H (i,c) vdW denotes the van der Waals interaction between the control atom in |1 c and Rydberg atom {|r ± i }.
As stated above, in our model, the effective spin-1/2 of the quantum simulator is encoded in ground states |g ± of the simulator atoms (see, for example, [36] and [43] for a specific choice). Following [36], an effective spin-spin interaction is obtained by admixing to these ground states with a laser a fine-structure split Rydberg state. We note that it is the combination of a finestructure resolved Rydberg manifold (i.e. the spin-orbit coupling) together with the van der Waals interaction (away from Förster resonances [44]) which provides the physical mechanism for the effective spin-spin interactions in the ground state manifold, and determines the spin models which can be realized in this setup (see Appendices A and B). To be specific, we assume below 87 Rb with Rydberg states |r ± = |nP 1/2 , m j = ±1/2 for the simulator atoms, and Rydberg state of the control atom |n S 1/2 , m j = 1/2 . In the basis {|r ± i ⊗ |r ± j } of the given Rydberg manifold, the interaction H (i,j) vdW is then represented by a 4 × 4 matrix [43], Here C 6 and C 6 are van der Waals interaction constants, I 4 is the identity matrix, and D 0 a 4 × 4 matrix with θ ij , φ ij angles representing the axis connecting the pair of atoms i, j with distance r ij . Explicit expressions for these quantities are provided in Appendix A 1. In a similar way, the van der Waals interaction between the simulator atom i and control atom c has the form  with C 6 the corresponding van der Waals coefficient. Here we made the assumption that C 6 is independent of the state |r ± (see Appendix A 2). Thus, conditional to the control atom to be in the Rydberg state |1 c , the Rydberg energies of the simulator atoms are shifted by E r± → E r± + C 6 /r 6 ic , which for the ring geometry of Fig. 1(b) provides an identical shift for all simulator atoms. In the spirit of the Rydberg blockade mechanism we assume this shift to be large. Thus, for a control atom in the ground state |0 c the dressing lasers ω ± are detuned from the Rydberg states |r ± i by ∆ ± , while the presence of a control atom in |1 c will detune these excited states by ∆ ± − C 6 /r 6 ic . For a control atom in |0 c the dressing lasers will thus induce interactions between the effective ground state spin-1/2 by admixing the Rydberg-Rydberg interactions H (i,j) vdW to the ground state manifold, while a control atom in |1 c detunes the Rydberg states, and thus effectively shuts off Rydberg dressing.

B. QND Hamiltonian and imperfections
The derivation of H QND proceeds now by a perturbative elimination of the Rydberg manifold of the simulator atoms. First, for the control atom in the ground state |0 c , and Ω c = 0, the dressing lasers for the simulator atoms are described by the detunings ∆ ± , as discussed above, and the relevant perturbation parameter for elimination of the Rydberg states is ξ ± = Ω ± /∆ ± 1. This yields an effective dynamics as spin-1/2 model for the dressed ground states |g ± i with effective spin-1/2 Hamiltonian denoted by H spin . The design of the desired spin Hamiltonians via dressing schemes by an appropriate choice of the atomic and laser configurations, and van der Waals interaction for a given geometry was discussed by [36]. In Appendix B we provide a derivation for the disordered 1D Heisenberg chain in the form (3). There we list explicit expressions for the J (x,y,z) , h etc. in Eqs. (B9-B12) of Appendix B as functions of the microscopic atomic parameters. Figure 2(a) shows the values of the spin-spin interaction constants J (η) ij as a function of the interatomic distance for a particular dressing scheme parameters.
Second, we repeat this derivation for the control atom in the excited state |1 c which amounts to replacing the detunings ∆ ± → ∆ ± − C 6 /r 6 ic in the expressions for the parameters of the Hamiltonian (B6) of Appendix B. This results in a spin Hamiltonian which we denote by H spin , with the same structure as H spin , but strongly suppressed couplings J J etc. After combining the two cases, we obtain the Hamiltonian describing the coupling between the simulator and the control qubit If H spin and H spin commute,H is equivalent to the QND Hamiltonian (1) with H spin → H spin − H spin . In the opposite case, H spin results in errors in the entangling gate U QND (t) , with the error rate characterized by a dimensionless parameter κ 3 = |H spin |/|H spin |, where | . . . | denotes the difference between the maximal and the minimal eigenvalues of the corresponding operator. Another source of error is related to spontaneous emission of Rydberg states of the control and simulator atoms, with error rates characterized by κ 1 ≡ γ d /|H spin | and κ 2 ≡ ξ 2 ± γ d L/|H spin |, respectively, where γ d and γ d are the corresponding spontaneous emission rates [see Fig. 2(b) and Appendix D for more details]. The maximal of κ i determines the coherence time t coh ∼ [|H spin | max{κ i }] −1 below which the errors in the gate operation due to spontaneous emission and H spin can be ignored.
Thus for t < t coh the quantum simulator coupled to the control qubit can be described by the QND Hamiltonian Eq. (1) which is written here in the rotating frame with respect to the laser of the control qubit, and we have set Ω c = 0 in H c (Ω c = 0 is needed only during the preparation stage of the control qubit and for the measurement readout).
In the above equation δ refers to a renormalized detuning of the control laser (see Appendix B).

III. MEASUREMENT OF THE SFF IN HAMILTONIAN SYSTEMS
For a generic quantum many-body system, the energy level statistics is an indicator which allows to distinguish between quantum ergodic and non-ergodic regimes [45][46][47][48]. Quantum ergodic (or chaotic) systems are characterized by eigenenergies statistics given by the Wigner-Dyson distribution [49] in a universality class of random matrices, and the eigenfunctions are delocalized over the configuration space. In contrast, non-ergodic quantum systems, such as integrable models, display Poisson statistics of energy levels and localized wave-functions. These considerations apply not only to systems with time-independent Hamiltonians and their energy spectrum, but also to timeperiodic Floquet systems [45] characterized by quasienergies.
An equivalent characterization of ergodic vs.non-ergodic dynamics is provided by the spectral form factor (SFF), where spectral features are displayed in the time domain, i.e. essentially as the Fourier transform of the two-point correlation function of the spectral density. For a manybody Hamiltonian H spin with eigenenergies E the SFF is defined as Here f (E) is a nonnegative normalized smooth filter function of width ∆E covering a band in the middle of the energy spectrum, where the density of states is flat, thus probing the properties of typical states of the many-body system, while eliminating contributions from spectral edges. The overline indicates a possible disorder average. The second line of (10) rewrites the SFF in terms of the microcanonical density operator ρ mc = f (E ) | | representing an initial density matrix. The SFF has been central in high-energy physics [29] and the recent discussion of the transition from many-body localization (MBL) to quantum chaos in a class of generic spin chains with disorder [30][31][32].

A. Properties of the SFF
To illustrate the generic behavior of SFF, we show in Fig. 3 a numerically calculated SFF for the disordered Heisenberg XXZ spin-1/2 chain of length L = 12 with Hamiltonian (3). A random local transverse field [last term in Eq. (3)] plays the role of disorder, with values h i ∈ [−W, W ] drawn from a uniform distribution and W characterizing the disorder strength. We assume periodic boundary conditions, and consider the sector with zero total spin projection on the z-axis, The choice of this model is based on the following considerations. First, it exhibits a transition from quantum chaos for weak disorder to many-body localization for strong disorder, and second, it belongs to the class of models, where the QND-gate and thus the SFF protocol of the following section can be implemented with the dressing scheme in Rydberg tweezer arrays (see Sec. II).
According to Fig. 3, for small times τ the SFF K(τ ) decays from its initial value K(0) = 1 due to dephasing on a time scale ∼ 1/∆E. Here ∆E ≈ (E max − E min )/6, i.e., 1/6 of the width of the energy spectrum. The signature of quantum chaos is the existence of a ramp at long times τ : Random Matrix Theory (RMT) for the Gaussian orthogonal ensemble (GOE), which is applicable to systems obeying time-reversal symmetry, predicts which is shown as dashed line in Fig. 3. In the formula above, τ H = 2π/δ E ∼ 2 L /L denotes the Heisenberg time associated with the mean level spacing δ E = E +1 − E in the middle of the spectrum. Furthermore, K ∞ = f 2 (E ), where the value of K ∞ is equivalently given by the inverse of the number N ∆E of the eigenstates in the energy interval ∆E, K ∞ ≈ N −1 ∆E . As shown in Fig. 3, the RMT prediction agrees well with simulations for a finite size chain in the limit of weak disorder (blue line) for times τ > τ Th , with τ Th the Thouless time defined as the onset of the ramp [31]. The flattening with increasing disorder strength (orange and green lines) is indicative of the crossover towards nonergodic (MBL) behavior.
The Hamiltonian (3) is time-reversal symmetric and thus its chaotic phase is described by the GOE. As discussed in Appendix C, this model can also be realized with complex J, so that the chaotic phase is described by a GUE. In systems with broken time-reversal symmetry, K(τ ) follows the RMT prediction for the Gaussian unitary ensemble (GUE), Our SFF protocol below exploits repeated QND measurements of the control qubit to both prepare the desired initial energy distribution f (E) in (10) [see Fig. 1(a)], as well as to read the SFF (10). Challenges faced in measuring the SFF, and addressed below, are decoherence times in quantum simulators, relative to times required to identify the 'ramp', and the number of measurements required to provide clear signatures of both the ergodic and non-ergodic regimes. Figure 3 plots simulated measurements for a finite measurement budget which compare favorably with the (exact) numerical results for the SFF (solid lines). The question to be addressed is whether, with given experimental resources, it is possible to see the main characteristics of the SFF.

B. Measurement Protocol via QND-Coupling to a Control Qubit
An experimental protocol to measure the SFF in the Hamiltonian case will require: (i) the ability to prepare an initial (microcanonical) ensemble of states ρ mc = f (E ) | | with a filter function f (E) of given width ∆E in the center of the spectrum; (ii) the ability to resolve for a given number of measurements, and thus signal-to-noise ratio the baseline K ∞ ≈ N −1 ∆E and the value of the SFF around the Thouless time; and finally (iii) the ability to observe for a given decoherence time (part of) the ramp τ Th < τ < τ H and possibly the long time behavior of K(τ ).

Preparation (verification) of a microcanonical ensemble
The first step of the protocol requires preparation of the microcanonical (MC) ensemble ρ mc = f (E ) | |. This can be achieved with a low resolution phase estimation algorithm (PEA), providing a probabilistic preparation of an energy band via repeated measurement of the control qubit. The PEA is based on the QND Hamiltonian (1), and requires a minimal number of measurements M of the control qubit to achieve a given measurement precision ∆E [20][21][22]. Figure 4 illustrates the idea of the preparation protocol. First, the many-body spin system initialized in a state ρ in and the control qubit prepared in the state |+ c are entangled by the QND interaction H QND (1) during a time t 0 . As a result the control qubit is rotated proportionally to the values of eigenenergies E of the spin Hamiltonian (3). The time t 0 is chosen to maximally spread the full spectrum E of the Hamiltonian H spin over the equator of the control qubit Bloch sphere. After that, the control qubit measurement in |± c basis shrinks the populated energy window of the spins state by a factor of ∼ 2. Repeating the cycle with increasing interaction times t m and postselecting "+" outcomes results in the state ρ mc with a narrow energy distribution. Figure 1(a) shows the quantum circuit for the full preparation procedure involving M = 3 control qubits entangled  with the simulator via Note, that the same result can be achieved by performing sequential entanglement and measurement cycles with a single control qubit. The interaction times are chosen as t m ≡ 2 m t 0 , m = 0, . . . , M − 1 where t 0 ≤ π/(|E − δ| max ). The renormalized detuning of the control laser δ introduced in Eq. (9) allows to tune the energy band of the state to be prepared. If we now select a run with all readouts "+ m ", m = 0, . . . , M − 1, then the initial state of the spin system is projected into the state ρ out with a narrow distribution in the energy eigenbasis | (see Appendix E) The success probability of the preparation is p mc ≡ p . The function (14) has a peak at x = 0 with a width ∼ 2 −M /t 0 , and, therefore, the conditional outcome state exhibits an exponential narrowing of the energy distribution/uncertainty around E ≈ δ. Note, however, that the variance of this distribution has the same scaling ∼ 2 −M /t 2 0 due to the presence of long tails. Following the above procedure, the microcanonical ensemble can be prepared by initializing the spin system in the infinite temperature state ρ in = ρ ∞ ∝ | | and postselecting the outcome state with the probability p mc ∼ 2 −M . As a result, the spin system will be probabilistically prepared in the microcanonical ensemble . The state has the mean energyĒ = Tr{H spin ρ mc } = δ and the bandwidth ∆E = Tr{(H spin − δ) 2 ρ mc } ∼ 2 −M/2 /t 0 . Note that the initial infinite temperature ensemble ρ ∞ can be sampled by random initialization of individual spins in up and down states where we additionally apply the constraint S z = 0 to probe the SFF in the zero-magnetization sector. An example of the resulting eigenstates probability distributions is shown in Fig. 5 for the Heisenberg chain (3) of L = 8 spins.
We remark that measuring the probability p mc of the successful MC ensemble preparation allows one to estimate (see Appendix H for details) the Heisenberg time τ H and the plateau value Here D is the Hilbert space dimension of the considered symmetry sector of the model. For the Heisenberg spin chain (3) with being the binomial coefficient (we assume L to be even). In the experiment the probability p mc is given by the ratio of the number of successful preparation attempts to the total number of experimental runs. The estimated values of τ H and K ∞ uniquely determine the behavior of the SFF in the chaotic regime assuming RMT (dashed lines in Figs. 3 and 6).
The low resolution PEA procedure can also be used to verify that a given state consists of a superposition of energy eigenstates in a narrow energy interval ∆E aroundĒ (MC ensemble): With δ =Ē, the appearance of M 2 log 2 (|E − δ| max /∆E) successive "+" readouts with probability close to 1 signals that the given state has a desired energy variance. Further, if the state ρ in contains a collection of excited states in a narrow energy interval created from some initial state by a time-dependent perturbation (see [33], for example), a MC state can be distilled with our procedure through several successive "+" readouts.

Protocol to measure the SFF
Following the preparation of the microcanonical ensemble ρ mc , we perform the evolution for a time τ with the QND-Hamiltonian (1), and finally measure the expectation values of σ x and σ y for the control qubit, as shown in the Fig. 5(a), . We obtain K(τ ) by repeating this sequence for different disorder realizations and averaging the result.
The SFF measurement scheme realizes a QND measurement meaning that the initial state of the spin system is not heated up or destroyed after the interaction with the control qubit. It is, therefore, possible to reuse the once prepared ρ mc for the measurement of K(τ ) with different times τ i (see Appendix G). The maximum number of recycling times N reuse is restricted by the coherence time t coh of the spin system as

Experimental challenges
An experimental realization of the SFF measurement faces two main challenges, which limit the achievable system sizes.
Time scales: -Propagation up to the Heisenberg time τ H , which grows exponentially with the system size L, is limited by the finite coherence time of the quantum simulator. Thus, observation of the behavior of K(τ ) at times τ ∼ τ H , requires coherence times t coh > τ H ∼ 2 L /(JL). We note, however, that the effects of interest such as the transition to chaotic dynamics at the Thouless time τ Th and the distinct behaviors of the SFF for quantum chaotic [K(τ ) ∼ τ ] and integrable systems [K(τ ) ∼ const.] take place at much shorter times to be compared with t coh (see Appendix D).
Signal magnitude: -The second challenge is the exponentially small values of the SFF ∼ 2 −L at the characteristic times. The threshold signal level which can be distinguished from the shot noise after averaging over N d realizations of disorder in the spin Hamiltonian and N measurements per one disorder realization is given by Thus, the total number of experimental runs per data point necessary to resolve the features of interest in K(τ ), is given by N run > 2 L √ N d /N reuse for the probabilistic preparation of the initial MC ensemble (see Appendix D).

Numerical simulation of the SFF measurement
To demonstrate the feasibility of the SFF measurement with our protocol, we perform numerical simulations of the measurement process in the disordered Heisenberg spin chain for various numbers L of spins. The simulation includes the probabilistic preparation of the initial state ρ mc and averaging over finite numbers of disorder realizations and experimental runs.
The results are presented in Figs. 3 and 6. According to Fig. 6, a small system of L = 8 spins with the measurement budget of N run 10 4 experimental runs per data point allows identification of the linear ramp in K(τ ). This indicates the chaotic behavior with level repulsion in the Heisenberg model with a weak disorder W = 2J (blue dots). In contrast, a strong disorder W = 10J results in the localized behavior of the system dynamics (orange dots) lacking correlations in the distribution of the eigenenergies. The colored lines and the corresponding shaded areas represent the numerical prediction for K(τ ) and the root-mean-square error of the simulated measurement, respectively. The black dashed line is the RMT prediction K GOE (τ ) Eq. (12).
The Thouless times τ Th for various disorder strengths can be probed in a larger system of L = 12 spins as shown in Fig. 3. It is evident from the figure that the time τ at which the data points approach the RMT prediction K GOE (τ ) given by the black dashed line grows with the increase of the disorder strength W , which is compatible with the expected behavior of the Thouless time τ Th . The root-mean-square error of the simulated measurement with N run 2 × 10 5 experimental runs per data point is shown with shaded areas around the solid lines indicating the numerical prediction for K(τ ).
In both Figs. 3 and 6 the horizontal and vertical gray lines mark the shot noise threshold K * and the coherence time t coh ∼ 10 2 J −1 in our Rydberg setup (see Appendix D), respectively. For the system of L = 8(12) spins we simulate MC ensemble preparation with M = 3(5) filtering steps, perform averaging over N d = 100 (20), and include recycling of the prepared MC ensemble for N reuse = 10 times. The parameters of the Heisenberg model (3) for both system sizes are ∆ = 1, J 2 = 0.1, It is important to stress, that in both cases of L = 8 and L = 12, the K GOE (τ ) curves (black dashed lines) are completely determined by the plateau value K ∞ for τ → ∞ and by the Heisenberg time τ H , both of which can be independently estimated as outlined in Sec. III B 1 above and discussed in Appendix H.

IV. MEASUREMENT OF THE SFF IN FLOQUET SYSTEMS
To study aspects of ergodicity, thermalization, and quantum chaotic dynamics, periodically driven or Floquet systems are particularly appealing for several reasons: First, due to the absence of energy conservation, generic Floquet systems can thermalize very rapidly and completely even for relatively small system sizes [50]. Further, since the density of quasienergies is generically flat, it is not necessary to unfold the spectra of Floquet systems to access their spectral statistics [45]. This relative simplicity of Floquet systems has led to intriguing recent results which, through explicit calculations of the SFF, establish an analytical connection between RMT and many-body quantum chaos in interacting periodically driven spin chains [30,[51][52][53]. As we discuss below, such Floquet spin models, and the experimental measurement of the SFF, can be realized naturally with Rydberg dressing schemes.

A. SFF of Floquet systems
In Floquet systems, we are interested in the statistics of eigenvalues of the unitary operator U (ϑ) which describes the evolution of the system during one driving period of duration ϑ: U (ϑ) | = e −iλ ϑ | with λ l ∈ [0, 2π/ϑ] being the quasienergy. Similar to the Hamiltonian case [see Eq. (10)], we define the SFF for Floquet systems as [51] where the integer "time" t is the number of the evolution periods, ρ in = f l | | is the initial quantum state, and the overline represents a possible average over disorder. Since generic interacting quantum chaotic Floquet systems exhibit a flat density of eigenstates in the thermodynamic limit, no spectral filtering or unfolding is necessary and one can use the infinite temperature ensemble ρ in = ρ ∞ = D −1 | | ∝ I as the initial state (here D is the dimension of the Hilbert space).
The behavior of the SFF in Floquet dynamics has the same characteristic features as in the Hamiltonian case. As an illustration, below we consider two Floquet systems. First, we present a Floquet model which exhibits a crossover between the circular orthogonal ensemble (COE) and the circular unitary one (CUE) as a function of the driving frequency ω = 2π/ϑ [54]. Then we consider kicked Ising models which demonstrate clear signatures of COE and CUE statistics even for a small system of L = 4 spins.
Crossover between COE and CUE statistics. -Here we consider an interesting example of a periodically driven system for which a random matrix class of the Floquet operator U (ϑ) differs from that of the time-dependent Hamiltonian H(t) which generates the dynamics of the system [54]. To be specific, we consider a piecewise constant Hamiltonian with H 1 during the first half of the driving period and H 2 during the second half. The evolution operator over one period is of the form where H 1 and H 2 are Heisenberg Hamiltonians with random magnetic fields h x,y,z i , which are normally distributed around zero with unit variance, The Heisenberg models can be realized in the quantum simulator based on dressed Rydberg atoms as discussed in Sec. II. For generic values of the driving period ϑ the spectral statistics of the Floquet operator U (ϑ) is described by the CUE. However, for finite system sizes and in the limit of high driving frequencies, the dynamics of the system is described by an effective time-independent Hamiltonian (H 1 + H 2 )/2 in which the random fields in the y direction cancel and time-reversal symmetry is restored. Therefore, in this limit, the spectral statistics of U (ϑ) belongs to the COE. We note that in this limit the density of quasienergies is determined by the effective Hamiltonian and is not flat. However, in the numerical examples below we find that effects due to a nonflat density of quasienergies are insignificant, and we take ρ in in Eq. (15) to be the infinite temperature ensemble also for high driving frequencies.
The SFF and the simulated measurement of the SFF (see Sec. IV B below) are shown in Fig. (7). Remarkably, even the relatively small Floquet system of L = 8 spins exhibits clear RMT behavior with a crossover between COE and CUE statistics upon changing the Floquet period ϑ. In particular, we observe COE statistics for short periods or high driving frequencies (blue dots) with ϑ < ϑ c ∼ 0.  for 0 < t < D where the Heisenberg time is set by the Hilbert space dimension D = 2 L , see Eq. (12). For long periods ϑ > ϑ c , the system is in the unitary class with CUE statistics (orange dots) resulting in K CUE (t) = t/D 2 for 0 < t < D (dashed line), see also Eq. (13). We note that for sufficiently large values of ϑ the initial decay due to dephasing lasts at most a few driving cycles. This is because for these values of ϑ the unitary U (ϑ) is not a sparse matrix, as it is for ϑ 1, but a dense one coupling practically all states in the Hilbert space with each other. As a result, chaotic behavior starts already at t = t Th ∼ 1.
Small Floquet systems. -Now we consider two kicked Ising models described by evolution operators where H x,y,z are the transverse Ising Hamiltonians with random magnetic fields h x,y,z i ∈ [−1, 1] The two models described by the Floquet operators U 2 (ϑ) and U 3 (ϑ) belong to the COE and CUE random matrix classes, respectively. Remarkably, the statistical distinction can already be seen in a small system of L = 4 spins as the simulated measurement of the SFF shows in Fig. 8.
In the quantum simulator based on Rydberg atoms the Ising models can be realized according to the general scheme presented in the Sec. II. These examples illustrate that, as in the case of the Hamiltonian systems, the SFF provides a sensitive tool for probing quantum chaotic behavior of Floquet systems. At the same time, the measurement of the SFF in Floquet systems faces the same challenges as are present in Hamiltonian systems, i.e., small values of the signal and exponentially long time scales. However, for experimental studies of many-body quantum chaos, Floquet systems can be beneficial because they typically exhibit pronounced RMT behavior even for comparatively small system sizes. Moreover, since spectral filtering is not required in Floquet systems, the preparation step described in Sec. III B 1 can be omitted.

B. Measurement Protocol
To generalize the measurement protocol for the SFF in Hamiltonian systems described in Sec. III to Floquet systems, we reformulate it in terms of the unitary operator U(t) which describes the coupled evolution of the system and the control qubit during t Floquet periods of duration ϑ, For simplicity, we focus here on Floquet systems with piecewise constant Hamiltonians H k for time periods τ k−1 < τ < τ k , where τ 0 = 0 and the time dependence is repeated periodically with period ϑ = k τ k . The corresponding Floquet operator reads U (ϑ) = k e −iτ k H k . The controlled evolution Eq. (19) is achieved by using the QND interaction Hamiltonian (1) between the spin system and the control qubit with H spin ≡ H k : The measurement protocol for the SFF starts with the initialization of the control qubit in the state |+ and the system spins in the infinite temperature state ρ ∞ . (In practice, it is sufficient to sample from the infinite temperature ensemble by preparing the system, e.g., in random product states.) Then, for a particular realization of disorder in the instantaneous Hamiltonians H k , we apply the controlled evolution U(t), and measure the expectation values of the operators σ x and σ y for the control qubit afterwards, Finally, the quantity | U (ϑ) t | 2 = σ x (t) 2 + σ y (t) 2 is averaged over disorders realizations resulting in the SFF K(t) for the Floquet system Eq. (15).
Experimental limitations of the measurement of the SFF in Floquet systems are analogous to the Hamiltonian case. As detailed in Sec. III B 2, these limitations originate from exponentially small values of the SFF in the presence of shot noise due to a finite number of experimental runs, as well as exponentially large values of D, the analog of the Heisenberg time, and finite coherence times in real quantum simulators. As an example, the effect of a finite number of measurements is illustrated in the Figs. 7 and 8 for the Floquet system described by Eq. (16) and Eqs. (17), (18), respectively.

V. CONCLUSIONS AND OUTLOOK
Recent experimental studies of ergodicity breaking in quantum many-body systems focus on the absence of thermalization of local observables in integrable [55][56][57] and many-body localized systems [58][59][60][61], and on the slow growth of entanglement [61][62][63]. While it is firmly established through theoretical work that key signatures of ergodic vs. nonergodic dynamics are carried by individual eigenstates of the Hamiltonian of a quantum many-body system and the statistics of the corresponding eigenvalues, accessing these signatures in experiments requires novel approaches beyond the standard paradigm of quantum simulation. As a first main result of the present work, we have developed a method that enables in-depth experimental studies of the level statistics of interacting quantum many-body systems through the measurement of the SSF, and we have discussed the feasibility of observing the SFF signatures of ergodic (RMT) vs. non-ergodic dynamics for both disordered 1D Heisenberg and Floquet spin models. We conclude that the key features of RMT in the SFF should be observable with presently available Rydberg experiments for system sizes of ten or more atoms. In a broader context, our method to measure the SFF is an example of a quantum protocol in which the state of a quantum simulator is prepared and monitored via the measurement of an auxiliary qubit, which is entangled with the quantum simulator by applying a QND gate. This is the second main result of our work: The implementation of a QND Hamiltonian H QND = H spin ⊗|0 c 0| with Rydberg tweezer arrays, which yields a QND gate as U QND (t) = e −iHQNDt . In the context of a many-body system engineered with a Rydberg tweezer platform, this QND Hamiltonian can be implemented for a broad class of spin models specified by a Hamiltonian H spin . While we consider 1D systems in form of a ring with the clock qubit in the center, our ideas also carry over to more complex 2D simulator geometries [1,4,36]. In addition, unique opportunities to combine quantum simulation with atomic clocks are offered by Alkaline Earth atoms [3,[64][65][66]].
An intriguing possibility opened up by the present study is the design and implementation of more general quantum protocols involving QND gates entangling the quantum simulator with a freely designable H spin with a set of control qubits. As noted before, this opens the door to run, for example, a quantum phase estimation algorithms on analog quantum simulators. In Ref. [67], a continuous readout of 'the energy' of a quantum many-body system was proposed as analog measurement of a homodyne current, with an implementation for a transverse Ising model with long range interaction. In contrast, quantum phase estimation based on the present QND gate with a freely designable H spin provides an essentially universal, digital quantum algorithm to achieve measurement and preparation of (a narrow band of) energy eigenstates in an analog quantum simulator setting. Moreover, quantum phase estimation can be utilized to compute the dynamical response functions of quantum many-body systems [68]. This opportunity is of particular interest in the context of NMR, where the ability to design H spin enables accessing the NMR spectra of molecules which are described by parametric Heisenberg models [69]. In this Appendix we provide details on the van der Waals interactions between nP 1/2 + nP 1/2 and nP 1/2 + n S 1/2 states as relevant for the model of Sec. II. Our discussion adapts and extends Ref. [36].
For any pair of atoms i and j, the dipole-dipole interaction Hamiltonian reads [71] where d (i) is the dipole operator of i-th atom and r ij is the relative distance between atoms. In second-order perturbation theory inV we obtain the effective van der Waals interaction [36,72] (we assume the absence of Förster resonances [44]) HereP is the projector onto the states of interest (nP 1/2 , nP 1/2 or nP 1/2 , n S 1/2 ), and Q βχ ≡ |β, χ β, χ| is the projector on manifolds that are populated only as virtual intermediate states, with δ βχ energy differences. We note that due to the perturbative nature of (A2), this expression is valid only beyond a certain critical radius r > r c (see below).

Interaction between simulator atoms
We first consider the van der Waals interaction between excited Rydberg states for the fine structure states |r α=± = nP 1/2 , m J = ±1/2 of simulator atoms. The relevant projector reads The possible interaction channels for the atoms in P 1/2 states are listed in Table I,  attributed to different scattering channels. (Explicit expressions involving Clebsch-Gordan and dipole matrix elements can be found in [36,72]), and corresponding plots for Rb atoms are shown in Fig. 9(a). The 4 × 4 matrix D 0 (θ, φ) referred to in Eq. (7) reads , and following the geometry of our setup in Fig. 1(b) we set θ = π/2. Thus we obtain for the interaction Hamiltonian the structure where W (i,j) The specific spin models which can be engineered via Rydberg dressing [36,37], i.e. by admixing the van der Waals interactions (A3) to the ground states by off-resonant laser light, is determined by the structure of the matrix elements (A4).

Interactions between simulator and control atoms
Here we consider the van der Waals interaction between the simulator and control atoms in the Rydberg states. We restrict the states of the control atom to the Zeeman manifold which includes the logical qubit state |1 c , i.e. |1 c = n S 1/2 , m j = 1/2 , and |1 c ≡ n S 1/2 , m j = −1/2 . As we discussed below, the unwanted coupling between |1 c and |1 c can be minimized by proper choice of the principal quantum number n [43]. For the states of simulator atoms we again consider The relevant scattering channels are listed in Table II. Similar to (A 1) above, we obtain for H (i,c) vdW in the basis (|r + |1 c , |r + |1 c , |r − |1 c , |r − |1 c ) an expression analogous to Eq. (7), where coefficients C 6 and C 6 are defined as above. Numerical results for C 6 , and the relative strength of the second (anisotropic) term in Eq. (A5) are shown in Figs. 10(a,b), respectively, for different principle quantum numbers n and n = n + ∆n. As can be seen from Fig. 10(b), the second term in V (i,c) vdW can be made much smaller than the first one by choosing n and ∆n properly. This makes the interaction essentially isotropic.
We note, however, that, even when the condition C 6 |C 6 | is satisfied, the second term in Eq. (A5) may cause unwanted transitions between the Rydberg states of the control atom, |1 c ↔ |1 c , which can bring the control atom out of the Hilbert space of interest. To suppress such transitions, we make them strongly off-resonant by imposing e.g. an external magnetic field. Thus we assume that the Hilbert space of the control atom can be represented by the two qubit states |0 c and |1 c , and we can describe the interaction between the simulator atom and the control qubit by the 'blockade' Hamiltonian (8).

Validity of the van der Waals Hamiltonian
Before proceeding, we emphasize that the validity of the van der Waals interaction Hamiltonian Eq. (A2) requires the conditions [71,73] to be satisfied: which sets the lower bound on the interatomic distance, r ij > r c , for i, j = 1, . . . , L, and r i,c > r c , where the critical distances r c and r c are defined as the distances when the left-hand-site in the above condition equals unity. The calculated dependencies of r c and r c on the main quantum numbers are shown in the inset of Fig. 9(a) and in Fig. 10(c), respectively. For the simulator atoms, however, direct diagonalization of V (i,j) dd , Eq. (A1), shows the presence of an avoided crossing at r ij ≈ 2.4µm, see Fig. 9(b). Therefore, to stay away from this nonperturbative situation, we take r c = 2.4µm and (using n = 60 and n = 71 as principal quantum numbers) for the shortest interatomic distances in our setup.

Appendix B: Derivation of the effective Hamiltonian
In this Appendix we derive the effective Hamiltonian of the simulator-qubit interaction which we show to be given by Eq. (9). We perform an adiabatic elimination of excited Rydberg states of simulator atoms |r ± following [36]. We first transform Hamiltonian H s [see Eqs. (4)-(6)] into rotating frame: where and The explicit expression of van der Waals interaction Hamiltonian according to Eq. (A3) is We now define the dressed ground |g α i and Rydberg |r α i states as the eigenstates of the Hamiltonian H where λ ± =[(∆ α /2) 2 + Ω 2 α ] 1/2 ± ∆ α /2 (we assume ∆ α > 0 and real Rabi frequencies Ω α ).
We now define the projectors onto the subspace of the dressed ground states |g α i Equivalently, we can write P g as a sum of projectors onto the eigenstates of the Hamilto- g is the projector onto the subspace formed by the eigenstate with the eigenenergy E (0) l . Due to the particular form of the blockade interaction Eq. (A5) we can perform the adiabatic elimination of excited states simultaneously for the two states of the qubit. We now derive the effective Hamiltonian of the simulator atoms using projectors, and write an effective Hamiltonian in the P g subspace up to two lowest orders in V as [74] H eff = P g H s P g + 1 2 After evaluating this expression analytically up to the 4-th order in ξ α ≡ Ω α /∆ α , we obtain [36] the effective Hamiltonian in general XYZ form: We now show that the model can be reduced to the Heisenberg XXZ model assumed in the main text, and provide the corresponding interaction coefficients J As a result of the adiabatic elimination procedure the qubit logical states |0 c and |1 c respectively acquire additional spin-independent shifts β ij and β ij (see below). Thus in the rotating frame of the control laser frequency, the Hamiltonian takes on the form Eq. (9) with an effective detuning δ = ij β ij − β ij + E cr − E cg .
In the case of a perfect blockade, where V and W refer to the van der Waals interaction derived in (A 1), and we assumed equal Rabi frequencies and detunings, such that ξ α = ξ. If the energy difference between the subspaces is much larger than the coupling |J (η) ij | |∆ B |, this term can be neglected, such that the resulting low-energy dynamics (with typical energies ∼ J) takes place in a the sector with a fixed S z = i σ z i . To conclude, we note that in Eq. (B6) one can also generate a position-dependent magnetic field, h z → h z i = h z + δh z i . This can be achieved by using position dependent Rabi frequencies, Ω is a uniformly distributed random number, we obtain δh z i ∼ J.

Appendix C: Implementation of complex flip-flop phases
We now provide a way to engineer complex flip-flop coefficients in our spin models. As a starting point, let us consider the dressing lasers with the Laguerre-Gaussian spatial mode profile corresponding to the angular momentum l ± (note that it is sufficient to have only one nonzero l α ). In this case, the dressing term in the Hamiltonian Eq. (4) becomes L k=1 α=± where φ k,± ≡ 2πkl ± /L. After performing the adiabatic elimination procedure as in Sec. B of this Appendix, we obtain the effective Hamiltonian of the form following form (we omit here the ∆ B -term) with the flip-flop amplitudes acquiring nonzero phases φ kj = i α (φ k,α − φ j,α ). We note, however, that the sum of the flip-flop phases in Eq. (C2) across the system is always equal to an integer number of 2π, having, therefore, no effect on spectral statistics. However, physically relevant phase of the flip-flop amplitudes can be generated using stroboscopic engineering. Let us consider the case when, during the stroboscopic period T = 2π/ω s , we use the plane-wave dressing lasers with the Rabi frequencies Ω α for the time 0 ≤ t 1 ≤ T , and the Laguerre-Gaussian dressing lasers with the Rabi frequencies Ω LG α for the time t 2 = T − t 1 . Then, for the stroboscopic frequency satisfying the condition h, J ω s ∆ ± , the effective Rabi frequencies at the atomic positions are where We see that, by varying the ratio t 1 /t 2 , one can generate values φ eff k,α ∈ [0, φ k,α ] for the effective phases. With the effective Rabi frequencies (C3), the resulting phases of the flip-flop terms are not multiples of 2π anymore, and the corresponding Hamiltonian (C2) belongs now to the unitary ensemble.

Appendix D: Experimental considerations for SFF
In this Appendix we discuss in more details the experimental challenges of the SFF protocol, which limit the achievable system sizes. We also present estimations of the coherence times and the available system sizes for our Rydberg tweezer implementation.
Time scales: -Propagation up to the Heisenberg time τ H ∼ 2 L /JL is limited by the finite coherence time of the quantum simulator. First, the preparation time (with M measurements) is t prep ≈ t 0 2 M ≤ π(JL) −1 2 M −1 (see Sec. III B 1). Therefore, the lower bound on the coherence time is t coh > t prep ≈ (JL) −1 2 M −1 . Observation of the behavior of K(τ ) at times τ ∼ τ H , requires coherence times t coh > τ H ∼ (JL) −1 2 L . In the case of a dominant individual single-spin decoherence with the dephasing rate γ d , the corresponding condition reads J/γ d > 2 L . This limits observation of K(τ ) at times τ ∼ τ H to moderate system sizes (see estimates for the Rydberg tweezer array below). We note, however, that, even when the Heisenberg time τ H is not accessible, the transition to chaotic dynamics at the Thouless time τ Th and the distinct behaviors of the SFF for quantum chaotic [K(τ ) ∼ τ ] and integrable systems [K(τ ) ∼ const] takes place at much shorter times to be compared with t coh . In general, systems with smaller L are characterized by shorter characteristic times, and require fewer experimental runs to resolve the key features (see discussion below). On the other hand, finite size effects tend to wash out the characteristic features of K(τ ) signaling the chaotic behavior.
Signal magnitude: -The scaling of typical SFF values can be estimated as K ∞ ∼ N −1 ∆E for the late time plateau value and as K(τ ≈ τ Th ) ∼ N −3/2 ∆E [75] for the SFF minimum preceding the Thouless time τ Th , here N ∆E 1 is the number of eigenstates in the initial MC ensemble ρ mc . Preparation of MC ensemble (see Sec. III B 1) using M filtering steps (control qubits) produces ρ mc with N ∆E ∼ 2 L−M eigenstates. Since the protocol has a success probability p mc ∼ 2 −M and the QND measurement scheme can exploit the prepared state several times (N reuse ) one can perform N measurements per data point in N (1) run N/(p mc N reuse ) experimental runs per one disorder realization. On the other hand, the threshold signal level which can be distinguished from the shot noise after averaging over N d realizations of disorder in the spin Hamiltonian and N measurements per one disorder is given by K * ∼ 1/(N √ N d ) (Appendix F). We find the necessary number of measurements using the condition K * ∼ K ∞ . Thus, the number of experimental runs N run = N d N (1) run per data point necessary to resolve the features of interest in K(τ ), is given by N run > 2 L √ N d /N reuse . The number of filtering steps M does not affect N run for the probabilistic preparation scheme, therefore, it is enough to use M ∼ 3 to eliminate the contribution of spectral edges in the SFF. However, it might be possible to use a semi-deterministic preparation scheme, e.g., by populating excited energy eigenstates in some energy interval by driving or quenching the system [33] followed by verification via M filtering steps with high probability (∼ 1) of success. In this case, the necessary number of experimental runs per data point can be improved N run > 2 L−M .
Rydberg tweezer implementation: -We conclude with a discussion of imperfections for the Rydberg tweezer implementation, in particular decoherence rates and the effectiveness of the Rydberg blockade. We also elaborate on geometrical limitations on interatomic distances which ensure the validity of our Rydberg dressing and put a constraint on the maximal size of the system.
We start with the discussion of the decoherence effects which limit the duration of the experiment. They originate from the finite lifetime of the Rydberg atomic states and from the error rate in the gate operation caused by non-perfect Rydberg blockade mechanism (the term with H spin in (9)). The former is characterized by two dimensionless parameters κ 1 ≡ γ d / |H spin | ∼ γ d /(JL) and κ 2 ≡ ξ 2 ± γ d L/ |H spin | ∼ ξ 2 ± γ d /J, for the control atom and the spin system, respectively. Here γ d and γ d are the spontaneous emission rates for the corresponding Rydberg states, and for the system atoms we take into account the collective enhancement (∼ L) of the spontaneous emission rate due to highly entangled nature of the manybody excited states. The factor ξ 2 ± in κ 2 represents the admixture of the Rydberg state as a result of the dressing. In a similar way, the error rate in the gate operation can be quantified by the parameter κ 3 ≡ H spin / |H spin | which can be estimated [see Eqs. (B10) and (B11)] as where R is the distance between the control atom and the system atoms and R b ≡ 6 |C 6 /∆ ± | is the Rydberg blockade radius with C 6 being the interaction constant [see Eq. (8)]. We note here the very high power in the above estimate such that κ 3 decreases very rapidly for R < R b , giving, for example, κ 3 ∼ 10 −3 for R/R b = 0.8. The largest of these parameters κ = max{κ i } sets the upper bound for the time during which the evolution is coherent and follows the ideal QND Hamiltonian (1), t coh ∼ (JLκ) −1 .
As an example let us consider the Rydberg states with n = 71 (control atom) n = 60 (system atom) and the following parameters of the dressing scheme: ∆ ± = −9 MHz and ξ ± = 0.2. The decoherence rates for the atoms are γ d ≈ 2π × 318Hz, γ d ≈ 2π × 406Hz. The parameters characterizing the spontaneous emission are then κ 1 ≈ 4.4 · 10 −4 , κ 2 = 1.6 · 10 −4 for L ∼ 10 atoms and the Rydberg blockade radius is R b ≈ 6.5µm (see Fig. 10). For κ 3 to be of the same order or smaller, one should have R ∼ R max = 0.75R b ≈ 5µm as an upper bound on the distance between the system and control atoms. In our system κ 3 is the largest decoherence parameter and it limits the coherence times to t coh ∼ 10 2 J −1 .
The constraint on the distance between the system atoms is related to the validity of our dressing scheme which is to say the validity of the perturbative approach: It has to be larger than some minimal value, R ij > r c = 2.4µm, see Appendix A 2. Within our ring setup [see Fig. 1(b)] in which the radius is given by R max and the separation between system atoms is limited by r c , simple geometrical considerations give the maximal number of system atoms which is with both constraints, L max π/arcsin[r c / (2R max )] = 12, the number used in our numerical simulations.
The full sequence of control qubit measurements with outcomes v ≡ {v m } defines the resulting unnormalized state of the spin system where M( v, δ) = M −1 m=0 M(v m , t m , δ). As a result, the unnormalized probability p for the eigenstates | with the eigenenergy E to appear in ρ out is where P m (±, x) = {1 ± cos(xt m )} /2.
If we now select a run with all readouts "+", v m = {+}, m = 0, . . . , M − 1, we obtain the output state with the narrow energy distribution p = P (+ 0 , . . . The success probability of the protocol is p mc ≡ p . The expressions (E3) and (E4) are used in Sec. III B 1.

Appendix F: Shot noise in SFF Measurements
Here we determine the shot noise in the SFF measurement. As described in Secs. III and IV, the measurement of the SFF involves the estimation of the expectation values σ x,y for the control qubit by averaging results of N measurements. The statistical properties of such averaging can be described by introducing N copies of the control qubit and the spin system and considering the fluctuations of the collective spin S x,y = N −1 N k=1 σ x,y k of the N control qubits. In particular, we consider each copy of the control qubit and spin system to be initialized in the states |+ = 1 √ 2 (|0 +|1 ) and ρ in , respectively, and entangled via a controlled unitary U(τ ) = U (τ ) ⊗ |0 0| + I ⊗ |1 1|. In the case Hamiltonian dynamics we have U(τ ) = e −iHQNDτ and, consequently, U (τ ) = e −iHspinτ . Since we are interested in the noise due to a finite number of measurements, in what follows we consider a specific and fixed realization of disorder in the Hamiltonian H spin .
The consecutive measurements of the control qubits yield an averaged outcome which is given by an eigenvalue m x,y of the observable S x,y . The statistical distribution of the outcomes m x,y is characterized by the corresponding moments of the collective spin observables with respect to the state ρ(τ ).
The second and fourth moments of S x read This result and similar expressions for S y allows us to express the SFF and its fluctuations (for a single disorder realization) as More precisely, the variance of the SFF estimation for N 1 reads var [K(τ )] = S x (τ ) 2 + S y (τ ) 2 2 − S 2 x (τ ) + S 2 y (τ ) 2 ≈ 4 N K(τ ) + 4 N 2 . The signal-to-noise-ratio is thus given by SNR ≡ K(τ )/ var [K(τ )]. The SNR grows linearly SNR ∼ K(τ )N/2 with N up to N ∼ 2(1 + √ 2)/K(τ ) where the SNR becomes 1. For a given number N of measurements and a fixed disorder, values of the SFF above the threshold of K (1) * ≡ 2(1 + √ 2)/N can thus be determined with an SNR that is larger than 1. The threshold value averaged over N d realization of disorder K * = K (1) * / √ N d is presented in Sec. III B 3 and shown as horizontal lines in Figs. 3 and 6. A further increase of the number of measurements results in a slower growth with SNR ∼ K(τ )N /2. ρ mc . Since the state ρ mc is diagonal in the energy basis it commutes with the QND Hamiltonian (1). Therefore, the state is not perturbed after averaging over measurement results for a certain time τ : ρ out = M(+, τ, δ)ρ mc M † (+, τ, δ) + M(−, τ, δ)ρ mc M † (−, τ, δ) = ρ mc .
Consequently, one can recycle the prepared microcanonical sate as long as the decoherence in the spin system is negligible i τ i t coh . Furthermore, the above applies to any initial state as only the diagonal part of its density matrix contribute to the SFF. This result is used in the Sec. III B 2 and III B 4.