Measurement Protocol for the Entanglement Spectrum of Cold Atoms

Entanglement, and, in particular the entanglement spectrum, plays a major role in characterizing many-body quantum systems. While there has been a surge of theoretical works on the subject, no experimental measurement has been performed to date because of the lack of an implementable measurement scheme. Here, we propose a measurement protocol to access the entanglement spectrum of many-body states in experiments with cold atoms in optical lattices. Our scheme effectively performs a Ramsey spectroscopy of the entanglement Hamiltonian and is based on the ability to produce several copies of the state under investigation together with the possibility to perform a global swap gate between two copies conditioned on the state of an auxiliary qubit. We show how the required conditional swap gate can be implemented with cold atoms, either by using Rydberg interactions or coupling the atoms to a cavity mode. We illustrate these ideas on a simple (extended) Bose-Hubbard model where such a measurement protocol reveals topological features of the Haldane phase.


I. INTRODUCTION
Nowadays, entanglement is a central concept in many branches of quantum physics ranging from quantum information [1] to condensed matter theory [2][3][4] and highenergy physics [5]. Given a many-body quantum state ρ, a fundamental quantity in characterizing the bipartite entanglement between two subsystems A and B is the corresponding reduced density operator of a partition, ρ A ≡ tr B {ρ} [6]. In particular, the entanglement spectrum (ES), the spectrum of ρ A , denoted as σ(ρ A ), has proven to be a powerful theoretical tool to analyze entanglement properties in a quantum information context [6][7][8] and has more recently also attracted much interest in condensed matter physics to characterize manybody quantum states. For example, as first pointed out by Haldane and Li [9], the ES can serve as fingerprint of topological order since it mimics the excitation spectrum of chiral edge modes [10][11][12][13][14][15]. In this context, the central spirit lies in the so-called bulk-edge correspondence of the ES [12,13], a manifestation of the holographic principle [3,5,11]. Moreover, the importance of the ES has been discussed in the context of tensor networks [11,15], quantum criticality [16], symmetrybreaking phases [11,17,18], and, most recently, manybody localization [19][20][21] and eigenstate thermalization [22]. At the same time, it has been pointed out in the literature [23] that the ES can contain nonuniversal features requiring caution in using it as a tool to locate phase * hannes.pichler@cfa.harvard.edu † gzhu123@umd.edu transitions, in particular, for symmetry-breaking phases. Given the enormous interest in the ES as an important theoretical concept and powerful numerical tool, an outstanding challenge, however, is the direct experimental measurement of the ES in many-body systems, where a full state tomography is prohibitively expensive or even impossible. At the same time, the rapid development of cold-atom experiments [24][25][26][27][28][29] has offered unique opportunities to access quantities related to entanglement in a many-body system [30][31][32][33]. A remarkable achievement in this context is the recent measurement of the second Renyi entropy (S 2 = − log Trρ 2 A ), i.e., the purity of the reduced density operator in cold-atom experiments in optical lattices [34,35] following the theoretical proposals in Refs. [30,31].
Motivated by these recent developments, we present a protocol on how the ES can be directly measured in a quantum simulator. In analogy to a many-body Ramsey interferometry [36][37][38][39], where the evolution of a manybody state is conditioned on an ancilla system, we develop a scheme where the conditional evolution of the many-body system is determined by a copy of the density operator, acting as the Hamiltonian [40]. This is achieved by a sequence of global swap operations between two copies of the system, controlled by the ancilla. The Fourier transform of the ancilla's Ramsey signal reveals the ES and the corresponding degeneracies. We discuss the physical implementation of this scheme in a cold-atom setup, relying on a combination of single-siteresolved addressing techniques in optical lattices [41,42] and dispersive interactions with an ancilla atom via the Rydberg blockade mechanism [43][44][45][46]. We illustrate our protocol with an extended Bose-Hubbard model contain-  1. (a) Circuit representation of the protocol to determine the spectrum of a density operator ρ. It consists of n stroboscopic steps (operations) that each involve an ancilla system and two copies of the state under investigation. (b) Each step can be constructed from three basic operations: tunnel coupling between neighboring copies (blue), a controlled (dispersive) phase shift for atoms in the lattice based on the state of the ancilla system (green), and rotations of the ancilla qubit (red). For a description of these operations, we refer to the main text [Eqs. (4), (5), and (7); see also ing the Haldane phase [47,48]. Our detection of the degeneracies of the "ground state" in ES and the entanglement gap serves as a signature [14] of the symmetryprotected topological phase [49].

II. PROTOCOL
To outline the key idea we consider two quantum systems, labeled 1 and 2, in the quantum states ρ 1 and ρ 2 , respectively. Here, the Hilbert spaces for system 1 and system 2 are assumed to be isomorphic. The basic building block of our protocol relies on the identity [40] (1) where S is the swap operator that interchanges the quantum state of systems 1 and 2: S|ψ 1 ⊗ |ψ 2 = |ψ 2 ⊗ |ψ 1 . The right-hand side of Eq. (1) describes a coherent evolution of system 1, which is not generated by a Hamiltonian but instead by the density operator of the second system, ρ 2 . Below, we show how one can reduce the required unitary U S ( ) = e −i S to a set of simpler operations that can be implemented in state-of-the-art experiments with cold atoms.
The central idea is now to repeatedly perform the operation U S ( ) and obtain a stroboscopic evolution of ρ 1 with ρ 2 according to Eq. (1), using a new copy of ρ 2 in each step. Therefore, after n steps, we obtain the map ρ 1 → E ρ2 (ρ 1 ) = e −in ρ2 ρ 1 e in ρ2 . Thus, ρ 2 takes the role of a Hamiltonian for system 1, which evolves "for a time" t n = n . Monitoring this evolution allows us to access the spectral properties of ρ 2 , e.g., via quantum phase estimation using an ancillary quantum computer [40]. Here we use a simple Ramsey technique instead and employ an ancillary system (with basis states |0 and |1 ) to control the application of U S ( ) [cf. Fig. 1(a)]: where |± = (|0 ± |1 )/ √ 2. For an ancilla initially prepared in the state |0 , the measurement of the operator Z = |0 0| − |1 1| after n such controlled stroboscopic steps gives Z n = 1 2 tr{e −i2tnρ2 ρ 1 + ρ 1 e i2tnρ2 }. This expression is valid for small time steps, such that t n 1. For the general expression, we refer to Appendix A. In particular, for ρ 1 = ρ 2 ≡ ρ, one gets The set of eigenvalues of ρ, {λ α }, can thus be extracted by a simple Fourier transform of the measurement signal for different n. Note that the choice of ρ 1 = ρ 2 = ρ is not fundamental, but it is a natural one for an experimental implementation and renders the protocol sensitive to the largest eigenvalues.
The basic building block of this protocol is the unitary U S ( ) = e −i S on the joint system of 1 and 2. We note that the recently implemented schemes to measure the Renyi entropy of cold atoms [34,35] require a measurement of the expectation value of the swap operator S. For a many-body system, S can be decomposed into a product of local swap operators, such that a measurement is possible by local operations only [31,32,50]. Here we aim for a more ambitious goal since we want to apply a unitary that is generated by the global swap operator. This is a nontrivial task since S, and, therefore, U S are highly non-local. Remarkably, the protocol outlined above can nevertheless be implemented with operations relying only on experimental tools already available in cold atom experiments as discussed below, such as controlled tunneling between neighboring lattice sites [34], local addressability of individual sites [41,42] and dispersive interactions based on the Rydberg blockade mechanism [43].

III. COLD-ATOM IMPLEMENTATION
While the protocol discussed above is completely general, we now present an implementation thereof in experiments with cold atoms in optical lattices. For concreteness, we consider bosons in a one-dimensional optical lattice described by a Bose-Hubbard model [51], but the protocol equally applies to two-dimensional systems. We are interested in the entanglement (spectrum) between different lattice sites of an arbitrary state ρ, i.e., the motional degrees of freedom of the system. A realization of the above Ramsey interferometer with cold atoms consists of the following steps: To realize the controlled phase shift, an off-resonant laser to a Rydberg state |s b is focused on the lattice sites in subsystem A of copy k. If the control atom is in the state |0 , the ac-Stark shift leads to an acquired phase (left); if it is in the Rydberg state |1 , the Rydberg blockade mechanism leads to a shift of the state |s b , suppressing the ac-Stark shift (right). (c) Illustration of the basic operations in step k to determine the spectrum of the reduced density operator of only the first three lattice sites. Note that the blockade radius rB exceeds the size of the lattice. The rotations of the control atom U correspond to Rabi pulses on a Rydberg transition. (d) Analogous setup for 2D systems. Several copies are created in different layers of a 3D lattice, with the ancillary atom trapped nearby. All steps in the protocol are the same as in the 1D case. We note that the subsystem A does not need to be contiguous.
(i) Preparation: n + 1 identical copies of the manybody state ρ are prepared. Such copies are naturally realized in experiments with cold atoms in optical lattices [34], e.g,. in parallel one-dimensional (1D) tubes. Once the states are produced, we freeze the motion along each lattice for the duration of the whole protocol by rapidly increasing the lattice depth and turning off interactions between the atoms, e.g., via Feshbach resonances [52]. In addition, an atom, whose internal states represent the ancilla qubit, is trapped in a separate lattice site, close to one of the tubes (Fig. 2). We initialize the ancilla atom in a stable state, representing |0 , and use a highly excited Rydberg state to represent |1 .
(ii) Stroboscopic steps: A single run of the protocol consists of n stroboscopic steps. The kth step (k = 1, . . . , n) involves only the ancilla atom and the atoms in the two adjacent lattices k and k + 1 (Fig. 1). We label the sites in each lattice by j = 1, . . . , M , and the bosonic annihilation operator in the kth copy by a j,k . In each single stroboscopic step, three types of operations are performed: (a) A tunnel coupling is induced between all the sites in lattice k and k + 1 [ Fig. 2(a)]. By lowering the potential barrier between the two lattices, e.g., using a superlattice, the atoms are allowed to tunnel between the two copies, as described by the Hamiltonian: [34]. After the time t BS = π/(4J) the potential is ramped up again, realizing the so-called "beam splitter", (b) The second type of process involves a controlled phase shift where a π phase is acquired by the atoms in sites j ∈ A of lattice k depending on the state of the ancilla atom.
To achieve this, we make use of the internal structure of the atoms in the two chains, which so far were treated as structureless bosons in a stable internal state |s a . In order to realize the controlled phase shift, we couple the atoms in subsystem A of lattice k to a Rydberg state |s b using a laser with Rabi frequency Ω, which is detuned by ∆ from the |s a ↔ |s b transition. The corresponding ac-Stark shift gives the required π-phase shift. However, if the control atom is in the Rydberg state |1 , the Rydberg blockade mechanism suppresses this process [see Fig. 2 On a more formal level, the corresponding Hamiltonian for this process reads Here, we denote the bosonic annihilation operator for an atom in the Rydberg state |s b on site j in the kth copy by b j,k . The second line of Eq. j,l ). In the fully blocked regime, for |V , one can adiabatically eliminate the state |s b and obtain, in second-order perturbation theory, the effective Hamiltonian, H eff = j∈A j,k a j,k . Note that the distance dependence of the Rydberg interactions (V (c) j ) as well as interactions between lattice atoms (V drop out in this regime. By applying H eff for a time t phase = π∆/Ω 2 , one obtains the unitary given in Eq. (5). Such Rydberg-blockade-based gates have already been studied and demonstrated experimentally [39,43,46]. We point out that the AC-Stark laser field has to be applied in a site-resolved way [41,42] only to the sites in A of copy k. (c) The third operation is a simple single-qubit rotation, as realized by addressing the control atom with a resonant laser of Rabi frequency Ω c for a time t = /Ω c [ Fig. 2(c)]. The rotation angle 1 determines the stroboscopic step size.
Combining these operations one obtains a single stroboscopic step by , as can be shown using S 2 = 1 [53,54]. Note that this operation differs from the one given in Eq. (2) by an additional swap, U step = (1⊗S)U step . This swap is convenient since it guarantees that in each stroboscopic step a new (unused) copy is coupled to the state that is propagated according to Eq. (1) while involving only neighboring copies (see Appendix A) .
(iii) Measurement of the ancilla atom in the Z basis, with outcomes ±1.
The average over many runs gives Z n , and the eigenvalues of the density operator can be extracted via a Fourier transform of Eq. (3) over the simulated time t n . Note that one can choose any subset of modes A by addressing the corresponding sites with the ac-Stark laser in step (ii-b). In particular, the addressability allows measuring the spectrum of any reduced state ρ A , as well as the spectrum of the entire state ρ.
We note that in systems with a conserved quantum number, such as the total number of atoms, it can be useful to access the ES in a quantum-number-resolved way [55]. In fact, our protocol can be easily adjusted to measure the ES in different number sectors [56], via a preselection of the copies by measuring the total number of atoms in subsystem B in all n copies before step (i). Beyond the conceptual asset of obtaining richer information about the entanglement structure, this also has the advantage of increasing the spectral resolution, as discussed below and pointed out in calculations of ES from quantum Monte Carlo simulations [57].
We further point out that our protocol is not limited to the implementation of the controlled phase shift using an ancilla atom and the Rydberg blockade mechanism. Alternatively one can place the atoms in an optical cavity [58,59] and use different photon number states as ancillary system. The different ac-Stark shift experienced by the atoms allows us to implement the controlled phase gate [60].
While in the above discussion we focussed on a onedimensional setting, it is easy to see that the protocol equally applies to other configurations such as twodimensional systems. Instead of preparing several copies in one-dimensional tubes, in step (i) the copies are, in this case, created in two-dimensional (2D) layers as indicated in Fig. 2(d). All other steps are identical to the ones described above. Experimentally, site-resolved addressing is more challenging in this 2D setting, as it involves also layer-resolved addressing. This could be achieved, for example, by combining quantum gas microscopes with individual addressing techniques using magnetic-field gradients [61]. Probing the ES in such 2D systems would allow us to diagnose, e.g., topological order in realizations of fractional quantum Hall states with cold atoms [62,63].

IV. ILLUSTRATION FOR THE HALDANE CHAIN
In the following we illustrate this protocol on the example of a Hubbard model. Here, we focus on an analysis of the simplest model in one dimension with symmetryprotected topological (SPT) order [49] and show that the protocol presented in Sec. III allows us to determine the largest eigenvalues of the ES, especially its gap and topological degeneracies. In particular, we consider an extended Bose-Hubbard model with nearest-neighbor interactions: where n j = a † j a j . As illustrated in Fig. 3(a), we consider a superlattice with a two-site unit cell with the hopping amplitude t j and the on-site potential ε j . The interact-ing part of the Hamiltonian contains the standard on-site interaction with strength U and nearest-neighbor terms with strength V j . Such extended Hubbard models have been realized recently in experiments with cold atoms with off-site interactions stemming from magnetic dipoledipole interactions [28]. The ground-state phase diagram of the Hamiltonian (8) can support a topologically nontrivial Haldane phase with a nonlocal string order parameter [64,65]. This can be most easily seen in the "hardcore" boson limit (U → ∞), where at most one boson can occupy a single lattice site. In this limit the problem can be mapped into a spin-1 2 chain by a † j → (−1) j (σ x j +iσ y j )/2 (with σ j the Pauli operators). For a proper choice of the Hubbard parameters (cf. Appendix C), the Hamiltonian (8) can be mapped into the alternating-bond Heisenberg chain [64][65][66] In the case of J, J > 0, the ground state of Eq. (9) displays alternating strong-weak antiferromagnetic (AF) bonds in the chain, as illustrated in Fig. 3(b) by solid and dashed bonds. In the limit of vanishing coupling on every second bond (J = 0), the spins dimerize in spin singlets [indicated by ellipses in Fig. 3(b)], while the two spins on the end of the chain are decoupled and hence form free edge states. This leads to a fourfold degeneracy of the ground state in an open chain, which survives in the thermodynamic limit also for finite J , where the system is analogous to the spin-1 AKLT model [48]. The Haldane phase in this regime falls into the category of SPT order and is protected by dihedral symmetries, time-reversal symmetry and inversion symmetry [64]. Its topological properties can be probed by a direct measurement of the nonlocal string order parameter [64,65]. Here we use the Haldane chain as an example to illustrate how our protocol reveals topological order in terms of the ES, which is a more generic detection tool, as it also applies to situations where no such string order parameter exists. The ES in the ground state shows a fourfold degeneracy if the chain is bipartitioned along a J bond [65] [cf. Fig. 3(c)]: If the size of subsystem A is larger than the edge-state correlation length, ρ A can effectively be decomposed into two parts, ρ A = ρ edge,A ⊗ ρ bulk,A , where ρ edge,A describes the state of the edge mode in subsystem A (i.e., on the left in Fig. 3) and ρ bulk,A the state of the bulk modes in subsystem A [49]. The edge mode ρ edge,A contributes a (trivial) factor of 2 to the degeneracy in the ES, while another factor of 2 stems from the entanglement in the bulk wave function ρ bulk,A . This degeneracy can be seen as a simple illustration of the bulk-edge correspondence in the ES [12,13]: The bipartition between A and B effectively frees edge spin(s) in A and forms an "entanglement edge mode" [Fig. 3(b)]. The gap in the ES, remaining finite in the thermodynamic limit, separates the physical and entanglement edge modes from bulk excitations and shrinks  A). We indicate the measurement uncertainty due to shot noise ∆Zn = (1 − Z 2 n )/N shot for N shot = 200 measurements per point. (b) Fouriertransformed Ramsey signal s(ν) = n Z n cos(2tnν), with the uncertainty ∆s ∼ n ∆Z 2 n . Each eigenvalue of the density operator gives rise to a Lorentzian whose width is determined by the size of the stroboscopic step. The product of the height and the width of the peak reveals that it stems from d = 4 eigenvalues at λ = 0.25. This is depicted in the inset by the circles, where the horizontal bars correspond to the width of the peak and indicate the spectral resolution. To observe the signature of the fourfold degeneracy around nmax ∼ 150 copies are necessary ( = 0.1). The position of the maximum, together with the height and the width of the resonance, allows one to determine the weight of the peak, i.e., the degeneracy (d) of the corresponding eigenvalue (see text). The inset shows the position of the resonance and its degeneracy d for the peak at λ ≈ 0.25 for the different values of = 0.1, 0.05, 0.025, 0.01. In all cases, one finds d ∼ 4, as expected for fourfold degenerate eigenvalues. The horizontal bar indicates the width of the peak, σα.
when the correlation length increases with J /J.
In Fig. 4, we illustrate how our protocol would measure the ES and detect both the degeneracy of the ground state and the degeneracy of the reduced system. We plot in Fig. 4(a) the evolution of Z n as a function of used copies, for different values of the stroboscopic step size , applied to a system described by a mixture of the four degenerate ground states of the Haldane chain [with J /J = 0.1, as in Fig. 3(c)], i.e., a thermal state at a temperature well below the excitation gap. As shown in Fig. 3(c), the ES is dominated by the fourfold-degenerate largest eigenvalue λ ≈ 1/4. The Ramsey signal Z n shows the corresponding characteristic oscillations as a function of the used copies n. The decay of the signal, which is not captured by Eq. (3), is due to the higher-order contributions for a finite step size (see Appendix A 3 for an analytical formula). In Fig. 4(b), we show the Fourier transform of this signal, s(ν) = n Z n cos(2t n ν). Because of the decay of the signal, it does not display sharp peaks at the eigenvalues of the density operator, but each eigenvalue λ α gives rise to a Lorentzian profile (with width σ α ∼ ), centered around it, s(ν) α λα 4 σα (ν−λα) 2 +σ 2 α . Since the eigenvalues appear as resonance frequencies, but also as the weight of the peak, one can determine the degeneracy d of an eigenvalue by the product of the height and the width of the corresponding peak i.e., the degeneracy d = 4s(λ α )σ α . This product is shown in the inset of Fig. 4(b), revealing the fourfold degeneracy of λ = 0.25. Note that, for a time step = 0.1, this degeneracy is already visible for a propagation time t nmax ∼ 15, i.e., n max ∼ 150 copies.

V. SPECTRAL RESOLUTION
In general, the spectral resolution δλ that is achievable with this protocol is determined by the total evolution time [δλ ∼ 1/( n max )]. For a given n max , the optimal step size scales as ∼ 1/ √ n max , maximizing the total integration time t nmax , while minimizing the decay of the Ramsey signal due to Trotter errors (see Appendix A 3). The corresponding spectral resolution scales as δλ ∼ 1/ √ n max . According to the Nyquist sampling theorem the signal Z n needs to be measured only at times t n with n m/(2 ) (m = 1, 2, . . . ), such that the total number of measurement points scales as ∼ 1/(2δλ). Note that beyond the Fourier analysis illustrated above, there exist powerful numerical techniques, based on Prony's algorithm that are tailored to extract frequencies of such damped exponentials [67]. Moreover, note that shot noise is not a severe limiting factor in our scheme. The corresponding uncertainty can be bounded by ∆Z n = (1 − Z 2 n )/N shot and ∆s n max /N shot respectively. Thus, the signal-to-noise ratio for the smallest resolvable peak (at λ ∼ ) can be estimated as ∆s/s ∼ 1/ √ N shot . This is in striking contrast to other approaches used to measure the ES based on measurements of Renyi entropies of different orders [31,57,68]. In these approaches, the spectrum is determined by the roots of the characteristic polynomial, whose coefficients can be expressed in terms of the moments of the density operator via Newton's identities. One of the main difficulties in such an approach is that root-finding algo-rithms can be extremely sensitive to noise in the coefficients of the characteristic polynomial, i.e., shot noise [69].

VI. CONCLUSION AND OUTLOOK
In this work, we showed how one can access the spectrum of an arbitrary (reduced) density operator of a many-body system by implementing a hardware protocol, instead of full quantum state tomography. While here we proposed a interferometric scheme that is within technological reach with upcoming cold-atom experiments, the scheme could be generalized to other quantum simulation platforms such as ion trap and circuit-QED architectures. Moreover, one can envisage it as the main building block of a full quantum principal component algorithm [40], where instead of a single Rydberg atom several ancilla qubits control the protocol, and the spectrum could be extracted via quantum phase estimation [1]. While here our analysis focused on bosonic atoms, we expect that analogous methods can also be employed for fermionic species [50]. Finally, we want to point out that the concept of density matrix exponentiation discussed here has applications beyond spectral estimation as explored recently in Refs. [70,71]. H.P., G.Z., and A.S. contributed equally to this work.

Appendix A: Protocol
Here, we outline the derivation of Eq. (3) of the main text and elaborate on the decomposition of a single stroboscopic step into processes that can be implemented in cold-atom experiments.

Ramsey interferometer
The unitary giving the Ramsey interferometer with n steps described in the main text can be compactly written in the form U Ramsey = n k=1 U (k,k+1) step , where we imply an ordering of the product defined as n k=1 A k = A n A n−1 · · · A 1 . Also note that we indicate the copies on which the operators act in the superscript; e.g. S (k,l) interchanges the quantum states in copy k and l and leaves the other copies invariant. As in the main text we suppressed this whenever there is no danger of confusion to simplify notation. Here, we want to prove Eq. (3), which states To this end we first note that one can write Since S (k,k+1) are unitary operators acting as the identity on the ancilla Hilbert space one can write the measurement result of the interferometer with n steps, Z n , as Using the relation one can calculate the trace by successively tracing out copy k = 2, 3, . . . , n (in this order). This gives Z n = 1 2 (tr {M n (ρ)} + c.c.), where we define the map M (X) = (cos(2 ) + i sin(2 )ρ)X + sin 2 ( )(X − ρtr {X}). One can separate the leading-order contribution in the 1 limit as tr {M n (ρ)} = tr {ρ(cos(2 ) + i sin(2 )ρ) n } + R n (ρ, ), where the remainder can be bounded as follows: Here we used |tr ρ(cos(2 ) + i sin(2 )ρ) k | ≤ 1 for k = 0, 1, 2, . . . . For → 0 and n such that t n = n = const, we get |R n (ρ, )| ≤ (e 2 tn − 1) → 0 for t n 1, such that To leading order, this is equal to Eq. (3) given in the main text.

Single stroboscopic step
Here we explicitly show that the gate sequence described in the main text indeed gives rise to the Ramsey interferometer leading to Eq. (3). To this end, we first note the following identity that can be checked with straightforward algebra (using S 2 = 1) U (k,l) c−swap , decomposing the operator (2) into a rotation of the ancilla, U , given in Eq. (7) and the controlled swap U (k,l) c−swap (also called Fredkin gate) which exchanges the quantum state in copies k and l based on the state of the ancilla qubit as (A4) Next we decompose the Fredkin gate into the three operations (i)-(iii) given in the main text. To this end we recall that for bosons the local swap operator, S is the gate given in Eq. (5), such that U (k,l) where we used [Ũ (k,l) BS , U ] = 0. Noting thatŨ 2 BS = S (up to an irrelevant phase) we obtain the gate sequence given in the main text Finally, note that formally the operatorŨ BS differs form the beam splitter given in Eq. (4) by local phase shifts. However, these local phase shifts are irrelevant in the presence of an atom-number super-selection rule, such that we can simply make the replacementŨ BS → U BS .

Exponential decay for finite
The expressions given in Appendix A 1 can be used to calculate the measurement outcome for any . However, it is instructive to consider a simpler circuit where the evolution with U S ( ) is only applied only to one arm of the interferometer, while in the other arm of the interferometer the state is left unchanged. One can straightforwardly adjust the protocol to such a setting by using a third internal level for the ancilla system. Such a modified scheme is more convenient for an analysis of effects of a finite since it allows for an analytical calculation of the measurement result to all orders in . This modified Ramsey interferometer will give a measurement result Z n = {E(n, )}, with E(n, ) = tr{U S (a,b) . Since the products of swaps can be written as a cyclic permutation, one can use identities given in Refs. [31,36] to express this in terms of moments of the density operator: One can easily note that this is a sum of damped exponentials, with σ α = − 1 2 log(cos 2 ( ) + λ 2 α sin 2 ( )) ≥ 0 and φ α = − arg cos( )−iλα sin( ) (cos 2 ( )+λ 2 α sin 2 ( )) 1/2 . A finite thus simply gives rise to a damped signal and a renormalization of the oscillation frequencies. Using Prony's analysis [67], one can extract the φ α 's and the σ α 's from the measurement record, which in turn determine λ α and its multiplicity from λ α = e −2σα −cos 2 ( ) sin 2 ( ) and λ α = tan(φα) tan( ) . While these identities giving damped exponentials hold only when the stroboscopic evolution is applied just to one arm of the interferometer, we found (numerically) similar behavior in the two-sided version presented in the main text, when is small. In particular, for the example studied in the main text the two versions are identical.

Appendix B: Experimental imperfections
In this section we analyze the robustness of the protocol with respect to errors and imperfect implementations of the individual steps of the protocol. We consider imperfect beam splitters due to (i) residual on-site interactions U (with U/J 1), i.e., j J(a † j,k a j,k+1 + a † j,k+1 a j,k ) + U 2 (a † j,k a † j,k a j,k a j,k + a † j,k+1 a † j,k+1 a j,k+1 a j,k+1 ) , and (ii) imperfect timing t BS → π/(4J)(1 + ξ BS ), where ξ BS 1 quantifies the deviation from the ideal case. Similarly, we model errors of the controlled phase shift by an error in the (controlled) phase, i.e., t phase → π∆/Ω 2 (1 + ξ phase ), and the error in the rotations of the ancilla qubit by → (1 + ξ rot ). In all cases, we model ξ BS , ξ phase , and ξ rot as random variables with zero mean and standard deviation σ BS , σ phase , and σ rot . The corresponding Ramsey signal for all four cases is shown in Fig. 5. The calculations are done for the same parameters as in Fig. 3; however, because of the large Hilbert space, we show only results for small systems of four lattice sites in each copy. For such small systems, the residual on-site interactions do not play a crucial role, as one can see in Fig. 5(a). This is not surprising as the number of particles in the subsystem is small. One expects that the deviation from the ideal case and the sensitivity to residual interactions increases with the number of particles, i.e., with the size of the subsystem. Among the other imperfections we find that our protocol is most sensitive to errors in the controlled phase gate. In the example shown in Fig. 5(c) the deviations from the ideal case are visible for fidelities f phase ≡ 1 − σ phase 0.99. While such fidelities have not yet been achieved with Rydberg atoms [46], theoretical calculations indicate that they are within reach [72]. Finally, random errors in the rotation angle of the ancilla qubit [ Fig. 5(d)] can be simply interpreted as dephasing and lead to a decay of the Ramsey signal. They thus have a similar effect as Trotter errors because of a finite time step (see Appendix A 3), and they eventually limit the maximal spectral resolution.
Appendix C: Extended-Bose-Hubbard model and alternating-bond XXZ spin chain We begin with the extended-Bose-Hubbard (EBH) model with alternating-strength nearest-neighbor hopping and interaction: Here t, t > 0 captures alternating nearest-neighbor hopping and V, V > 0 captures alternating nearest-neighbor repulsive interactions, while U and ε, ε represent on-site repulsion and on-site potential energy respectively. In the infinite-U limit, i.e. "hard-core" boson regime, the model we consider here can be mapped into a spin-1 2 chain, by setting a † j → (−1) j σ + j and a † j a j → 1 2 (σ z j + 1). The definition of spin operators with an alternating sign (equivalent to a gauge transformation) is to remove the negative signs in front of t and t . At half filling (adjusting ε = V and ε = V ), one can get the alternating-bond XXZ model [64][65][66] described by the following Hamiltonian: where we have defined J = 1 2 t and J = 1 2 t . Because of the open boundary conditions one needs an extra on-site potential in the first and last sites (j = 1 and j = N ) of δε = V /2 for the mapping to be exact. The Haldane phase and fourfold degeneracy in the "ground state" of the ES exist in this model and can survive even in the limit of zero zz coupling, i.e. V, V = 0 [65]. For simplicity, in the main text, we have chosen the off-site interac-tions to be V = 4J and V = 4J such that we realize an isotropic alternating-bond Heisenberg model. However, we note that such fine-tuning is not essential for the observation of the degeneracy in the ES and its relation to topological order.
The ground states for J < 0 and 0 < J < J are in the same Haldane phase and can be adiabatically connected. At the quantum critical point J = J = 1, one recovers the Heisenberg model which leads to a gapless phase. For J > J, one enters to another Haldane phase with twofold entanglement degeneracy if we cut on the J bond and a fourfold entanglement degeneracy if we cut on the J bond.
In the XXZ spin model the total magnetization along S z = j σ z j is a good quantum number, [H eff , S z ] = 0, corresponding to the conserved number of atoms in the extended Hubbard model, N tot = j n j . Therefore, in the absence of external magnetic field, the fourfold degeneracy can be divided into different S z (N tot ) sectors. In the simple limit J = 0, the edge modes just involve the two spins in the end of the chain and are composed of the four basis states | ↑ L | ↑ R , | ↑ L | ↓ R , | ↓ L | ↑ R and | ↓ L | ↓ R . Here L/R refer to the left/right edge mode at site j = 1 and j = N , respectively. When mapped to bosons for the cold-atom realization, the four basis states become |1 L |1 R , |1 L |0 R , |0 L |1 R and |0 L |0 R . The four states can obviously be divided into three different number sectors, with total boson number on the edges being n = 0, 1, and 2. While the n = 0 and n = 2 sectors are not degenerate, the n = 1 sector is twofold degenerate. The same counting applies to the case when J = 0 and the edge mode becomes extended. Therefore, in the cold-atom experiments, if the total number of particles is fixed to N/2, the ground-state degeneracy is actually twofold at exactly half-filling. Nevertheless, in this case the entanglement degeneracy is still fourfold because of the contribution from both the physical and entanglement edge modes.
Since the SPT phase is short-range entangled, a disentangling into product states is always possible with a local unitary transformation or, equivalently, a finitedepth quantum circuit [49]. Therefore, as long as the system size is much larger than the correlation length, we are always allowed to decouple the edge and bulk states, i.e. ρ A = ρ edge,A ⊗ ρ bulk,A , after a local unitary transformation.