Quantum computation of dynamical quantum phase transitions and entanglement tomography in a lattice gauge theory

Strongly-coupled gauge theories far from equilibrium may exhibit unique features that could illuminate the physics of the early universe and of hadron and ion colliders. Studying real-time phenomena has proven challenging with classical-simulation methods, but is a natural application of quantum simulation. To demonstrate this prospect, we quantum compute non-equal time correlation functions and perform entanglement tomography of non-equilibrium states of a simple lattice gauge theory, the Schwinger model, using a trapped-ion quantum computer by IonQ Inc. As an ideal target for near-term devices, a recently-predicted [Zache et al., Phys. Rev. Lett. 122, 050403 (2019)] dynamical quantum phase transition in this model is studied by preparing, quenching, and tracking the subsequent non-equilibrium dynamics in three ways: i) overlap echos signaling dynamical transitions, ii) non-equal time correlation functions with an underlying topological nature, and iii) the entanglement structure of non-equilibrium states, including entanglement Hamiltonians. These results constitute the first observation of a dynamical quantum phase transition in a lattice gauge theory on a quantum computer, and are a first step toward investigating topological phenomena in nuclear and high-energy physics using quantum technologies.


I. INTRODUCTION
Understanding how strongly-coupled quantum manybody systems behave far from equilibrium is a common goal across many physics disciplines. In nuclearand high-energy physics, evolution of matter and the mechanisms for its thermalization and hydrodynamiza--This manuscript has been authored by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725 with the U.S. Department of Energy (DOE). The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paidup, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for the United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE's Public Access Plan.
In particular, focusing on far-from-equilibrium dynamics after a quantum quench, we investigate dynamical quantum phase transitions (DQPTs), first proposed and observed in a range of quantum many-body systems in Refs. [140][141][142][143][144][145][146][147][148], and later predicted to occur in gauge theories as well [149], see also Refs. [150][151][152][153]. Importantly, this phenomenon is shown to be an ideal candidate for exploration with existing analog and digital devices, since it persists in small systems and on short time scales [149,154,155]. DQPTs are manifest in non-equal time correlation functions and non-equal time wavefunction overlaps (Loschmidt echos), and in the Schwinger model they are related to an underlying topological transition [149]. In fact, from a phenomenological standpoint, there may be interesting connections between potential non-equilibrium dynamics that could change the value of the topological θ angle in the early universe and the strong CP problem, a possibility that quantum simulations of real-time phenomena in simple models, and eventually in QCD, may shed light on.
This manuscript is organized as follows: Sec. II contains i) a brief overview of the topological DQPT exhibited by the Schwinger model, our strategy for ii) simulating the system on a quantum computer, and iii) identifying the topological phase transition by analyzing dynamical non-equal time correlation-function holonomy. In Sec. III, entanglement tomography is employed to build a classical representation of the systems' entanglement structure for non-equilibrium states. Our results are summarized in Sec. IV. The manuscript is supplemented by several Appendices. Appendices A and B include analytical details and quantum algorithms for real-time evolution of the Schwinger model, including the scheme to obtain non-equal time observables. Appendix C contains experimental details of the IonQ device used in this study. Appendices D and F provide a description of a symmetry-based error-mitigation scheme employed in this work, and a discussion of device imperfections. Appendix G contains more details on the entanglementtomography results. All circuits, raw data and the data shown in figures can be found at https://gitlab.com/ Niklas-Mueller1988/dqpt_2210.03089.
on a one-dimensional (spatial) lattice with N sites and periodic boundary conditions (PBCs), with mass m, electric coupling e, and lattice spacing a. ψ † n and ψ n denote creation and annihilation operators for the (staggered) fermions, respectively. U n is the link and E n the electric-field operator, satisfying the commutation relation [E n , U m ] = δ nm U n . The Hamiltonian commutes with Gauss's law operator at each site, G n ≡ E n − E n−1 − Q n , where is the staggered fermion charge. The gauge-invariant physical Hilbert space contains states that satisfy G n |ψ⟩ phys = 0.

Non-equal Time Observables Entanglement Tomography
Non-equal time correlator g q (t)

State Preparation Time Evolution Quench
Hamiltonian evolution  (5)) and NECFs gq(t) (Eq. (7)), include a symmetry-based error-mitigation scheme. (c) entanglement tomography scheme to extract Rényi entropies, fidelities, as well as the reduced density matrix ρA(t) from an entanglement Hamiltonian ansatz that is constrained by a classical optimization based on a number of random measurements.
Adding a topological θ term eθ 2π n E n to Eq. (1), and upon performing a chiral transformation to absorb the θ parameter into a (complexified) fermion mass term, yields the Hamiltonian [183][184][185][186][187][188] H(m, θ) = This model is considered a prototype model for CP violation in QCD [189][190][191]. With this motivation, a quench of the θ parameter, mimicking a rapidly changing background axion field [192][193][194], was considered in Ref. [149], and a novel non-equilibrium topological transition was discovered in the far-from-equilibrium response of the system. The transition was also identified as a manifestation of the more general phenomenon of DQPTs [140-144, 146, 147, 150], and recognized as an ideal target for quantum simulators and computers because of short time scales involved and potential for its realization in small systems. This phenomenon will be explored in this work using a digital trapped-ion quantum computer.
In particular, we employ the maximum possible quench of the θ parameter of the Schwinger model, ∆θ = π, corresponding to changing the sign of the mass term. This is achieved by preparing the ground state |Ψ(t = 0)⟩ ≡ |GS(m)⟩ of H(m) ≡H(m, θ = 0) and then time evolving according to H(−m) ≡H(m, θ = π), Two different non-equal time observables will be computed. The first is the Loschmidt echo, the overlap of the time-evolved state with the initial state, from which one can define an intensive rate function, The non-analyticities of Eq. (6) correspond to DQPTs and can be extracted on lattices as small as four to eight sites with small finite-volume effects in the e/|m| ≲ 1 regime [149]. The second quantity is a set of non-equal time correlation functions (NECFs), defined in the staggered lattice formulation as where q ∈ [− N 4 , N 4 − 1] and with ⟨. . . ⟩ ≡ ⟨GS(m)| . . . |GS(m)⟩ and U n,m (t) ≡ n−1 k=m U k (t) (U n,n = 1). From g q (t), an integer-valued topological order parameter can be extracted [195] where Here, z ≡ (q, t ′ ), g z ≡ g q (t ′ ), andg z ≡ g z /|g z |. Concretely, C ± (t) runs clockwise (counter-clockwise) in the positive (negative) half of wavenumbers q in the (q, t ′ )plane, i.e., C + (t) : and similarly for C − (t). Note that t ′ is continuous and q is discrete, hence integral and derivative become a sum and finite difference along the q-sections of C ± (t). Equation (9), which is valid at arbitrary coupling e, changes by an integer whenever the system undergoes a dynamical quantum phase transition. We refer the reader to Ref. [149] for more details, and here we focus on a quantum-computational representation of the topological phenomenon. The gauge degrees of freedom can be integrated out using Gauss's law in a way that maintains translation invariance with PBCs, see Ref. [196] for details. A gaugefield zero mode (average electric field) remains untreated in this case, but since our goal is to demonstrate the quantum computation of topological NECFs and utilizing entanglement tomography on an existing device, instead of engineering gauge-field digitizations, we leave the treatment of this term to future studies [197].
Explicitly, in the purely fermionic form, U n is set to one in Eqs. (1) and (8), and the electric-field term in Eq. (1) is replaced by a (translation-invariant) long-range fermionic density-density interaction [149,196] where d nm = min(|n − m|, N − |n − m|), and The resulting Hamiltonian reads with H I defined above and being the non-interacting fermionic part of the Hamiltonian.
The fermionic degrees of freedom can be represented in two different bases. Firstly, they can be represented by fermionic Fock states in position space which are given by spin states, identifying an occupied (unoccupied) state with |η n = 1⟩ ≡ |1 n ⟩ = |↑⟩ (|η n = 0⟩ ≡ |0 n ⟩ = |↓⟩), where n ∈ [0, N − 1] labels a lattice site. The operator-ordering convention is adopted from left to right in ascending order, i.e., Here, Fock states are defined with the same convention as in position space, with |↑⟩ (|↓⟩) identified as the occupied (unoccupied) state, and with modes ascending in q, Note that the local momentum (q-) Hilbert space has four states, In this manuscript, fermionic parity is explicitly realized within the unitary operations; our mode ordering is exactly equivalent to a Jordan-Wigner transformation starting at site n = 0 (momentum q = − N 4 ), i.e., ψ . Position-and momentum-space computational bases are related by composite fermionic Bogoliubov and staggered Fourier transforms [198][199][200][201][202] (n ∈ [0, N/2 − 1] labeling a supercell with n = 2n for even sites, n = 2n + 1 for odd sites), where Ψ n ≡ (ψ even n , ψ odd n ) T and Ψ q ≡ (ψ even q , ψ odd q ) T . Here, the even (odd) superscript on the position-space fields denote that they are acting on even(odd)-indexed sites n, while even (odd) superscripts on momentumspace fields indicate that the fields are obtained from their even (odd) position-space counterparts. The Bogoliubov transformation transforms creation and annihi-lation operators as follows with To uncover the DQPT, we study the non-equilibrium dynamics after a quench of the topological angle by ∆θ = π, i.e., from +m to −m, by quantum computing nonequal time observables, focusing first on the Loschmidt echo L(t) and the rate function Γ(t), defined in Eqs. (5)(6). Our general approach is summarized in Fig. 1 (a) and the schemes to measure non-equal time observables are shown in (b).
To do so, first the (non-interacting) ground state of q ⟩ are zero-fermion number eigenstates of H 0 (+m), is prepared in a computational-basis state in the (prequenched) momentum-space representation (a bar distinguishes them from eigenstates of the postquench H 0 (−m)). Then a Ramsey interferometry scheme is applied, depicted in Fig. 1(b) and discussed in greater detail in App. B 2, with an ancilla in Hadamard superposition representing the two interferometric paths of the qubits encoding the fermions. We begin with the non-interacting case, e = 0, where the initial state is time evolved in momentum space with no Trotter error. Since the system is prepared in the ground state of H 0 (+m) with positive mass, +m (Eq. (15)), but is evolved with H 0 (−m) with negative mass, −m, and the dispersion ω q is independent of the sign of the mass, H 0 (+m) and H 0 (−m) are formally identical in momentum space and the quench is not realized by simply changing a parameter in H 0 . Instead, a change of basis, from H 0 (+m) to H 0 (−m), is performed, see Appendix A for details. At e/|m| > 0, a Trotter scheme is applied, where the time-evolution operator is split into a free part (H 0 in Eq. (15)) realized in momentum space, and an interacting part (H I in Eq. (11)) realized in position space, where δt ≡ t/N T , with N T being the number of Trotter steps. V ≡ F B contains the fermionic Fourier, F , and Bogoliubov, B, transforms defined in Eqs. (16)(17)(18).
In their respective bases, the time-evolution operators e −iH I δt and e −iH0(−m)δt are diagonal and their control by the ancilla is easy to implement, see Appendix B for details. The basis transformations V and V † need not be controlled. The quench in the mass parameter is still performed via a basis transformation [203]  produced by the quantum computation. Errorbars in (a) are from statistical uncertainty assuming a binomial distribution for the outcome of single-qubit measurements, see Appendix E. The same uncertainty is propagated to (b) where, because of the logarithm in Eq. (6), deviations are largest in the peak region. An error-mitigation scheme is used in this work, relying on the measurement of the system qubits in addition to the ancilla qubit. Because of the particle-number conservation in the model, events where the measured system bit sequence violates this constraint can be identified and eliminated at the cost of reduced statistics for a fixed number of shots. The bottom panel of Fig. 2 shows the fraction of physical (symmetry-preserving) results (red) versus all results (gray), with a symmetry-violating device error occurring in about 25% of the events. Machine errors are not fully eliminated by this procedure, as e.g., errors that conserve particle number are not mitigated. Details, including a comparison of the unmitigated and error-mitigated data can be found in Appendix D.
Both the Loschmidt echo L(t) and the rate function Γ(t) for N = 8 lattice sites are shown in Fig. 3. To resolve the DQPT, significantly more shots, up to n shots = 16, 000, are required in the peak region at t|m| ≈ 1.1 (gray dashed line). The quantum hardware performs to the extent that a feature is discernible, reproducing, with some deviation, the simulator results after error mitigation. Errorbars are from the statistical shot-noise uncertainty effecting Γ(t) mostly in the peak region.
Because it is more difficult to prepare the ground state of the interacting theory, e/|m| > 0, a double quench is considered: First, the ground state of the non-interacting theory with H(m, e = 0) is transformed to that with H(−m, e = 0), then time evolution is performed with H(−m, e), at finite e/|m|. Using the Trotter scheme, Eq. (20) and given basis transformations results in a significantly greater circuit complexity compared to the non-interacting case. Figure 4 investigates the effect of Trotterization. The rate function Γ(t), computed exactly, as a function of time t|m| and coupling e/|m| is shown for (a) N = 4, |m|a = 0.9 and (b) N = 8, |m|a = 0.8. The position of the DQPT at t = t c is further shown, indicated by the change of the topological order parameter ν(t) in an integer step from ν = 0 to ν = 2, for the exact time evolution (black lines), with 1 (orange dotted lines) or 2 (purple dashed lines) Trotter steps. While agreement is expected at e = 0, where the Trotter error vanishes, the transition can be quantitatively reproduced even at large couplings e/|m| ≳ 1 with only one Trotter step. Figure 5 summarizes the results of quantum computing (a) L(t) and (b) Γ(t) at e/|m| = 1, N = 4, and with N T = 1 Trotter step. IonQ results (red symbols) are noisier when compared to ideal-simulator results (blue symbols), due to an approximately four-fold larger gate count relative to e = 0 [204].
Symmetry-based error-mitigated results are shown in Fig. 5, and the bottom panel of figure shows the fraction of physical (red bars) compared to all measurements (gray bars), with about 70% of events containing an identifiable symmetry-violating error. Unmitigated errors are significant, and the position of the DQPT is not statistically confirmed with the presented data set. Coherence and error models, which could potentially robustly resolve the DQPT and optimize this calculation, are discussed in Appendix F. We show statistical shotnoise uncertainty, noting a severe degradation of the machine results (red) because of the reduced statistics due to symmetry-based postselection.

B. Quantum computing topological transitions
A time-dependent topological index ν(t), as defined in Eq. (9), is associated with the DQPT in this theory, as discussed in Sec. II and in Ref. [149]. Errorbars are from the statistical shot-noise uncertainty which is nearly negligible. The topological index can be obtained by quantum computing the non-equal time correlation functions g q (t) in Eq. (7) using Eqs. (9-10). This is done by applying a Ramsey interferometry scheme using an ancilla in Hadamard superposition, Fig. 1(b), involving (uncontrolled) forward time evolution, followed by the ancilla controlling an annihilation operator in position space, ψ n , followed by (uncontrolled) backward time evolution, and the control of a creation operator in position space, ψ † n , see Appendix B for details.
Because of the significant gate cost at e ̸ = 0 resulting in reduced performance, we focus on e = 0 where g q (t) can be computed separately for each q, hence a lower gate cost. Explicitly, where the unitarity of U q in Eq. (18) is used. Here, ⟨·⟩ ≡ ⟨GS(m)| · |GS(m)⟩ refers to the expectation values in ground state of H 0 with mass m, while a q and b −q are Fock operators with respect to the eigenstates of H 0 (−m) in momentum space. We make use of a q-local Jordan-Wigner transformation to map fermion creation and annihilation operators onto spin raising and lowering operators, a † q = σ + q,a , and b † −q = −σ z q,a σ + q,b . Because spin raising and lowering operators are not unitary, they cannot directly be realized in the circuit and are decomposed into unitaries σ x = σ + + σ − and σ y = i(−σ + + σ − ), with eight combinations ⟨σ x/y q,a (t)σ x/y q,a (0)⟩ and ⟨σ x/y q,b (t)σ x/y q,b (0)⟩ which are computed separately [205]. The circuits use three qubits to compute real and imaginary parts of Eq. (21) (for both target and neighboring "spectator" qubits (which occur typically due to spillover of control signals) as well as a relative increase in ion heating, motional noise, mode frequency drift, and axial motion due to parallel operations. Nonetheless, since the circuits here are modest six-qubit circuits with relatively shallow gate depth, we expect no significant differences compared to separate executions of three-qubit circuits, hence adopting this computationalcost-saving option.
The results of quantum computing the NECF, g q (t), for (a) N = 4 and (b) N = 8 lattice sites and zero coupling e = 0 are displayed in Fig. 6. Exact results for both real (solid lines) and imaginary (dashed lines) parts of g q (t) are shown in the plots, along with the simulator (blue symbols) and IonQ results (red symbols), the latter computed with n shots = 500 shots. Exact results are reproduced well at N = 4 and N = 8, for all q-modes. A symmetry-based error-mitigation scheme is applied, where results that violate a conservation law are dropped, see the bottom panels of Fig. 6(a) and (b) [206]. Because of the small system size and the relatively small unphysical part of the 2-qubit Hilbert space, the noisemitigation scheme has little effect. Errorbars are from the statistical shot-noise uncertainty.
Furthermore, the time-dependent topological index ν(t), extracted from g q (t) via Eqs. (9)(10), is plotted in and similarly S B defined in Eq. (22)), for N = 4, |m|a = 0.9, e = 0, nCUE = 25, and n shots = 1000, including simulator (blue symbols) and IonQ (red symbols) results. The middle panel depicts fidelity F(t) (Eq. (23)). The bottom panel shows the Rényi entropy of the full system, S A+B (relative to the environment). Blue horizontal lines in the middle and bottom panels indicate ideal results, and a horizontal red line indicates zero fidelity or maximal entropy (ρ(t) = I/2 N ), respectively.
pute Eq. (10), the data in Fig. 6 is linearly interpolated in time. The topological index ν(t) is successfully reconstructed from the IonQ data (red symbols) for (a) N = 4 and (b) N = 8 [207] and reliably indicates the position of the DQPT. Because ν(t) is topological, essentially measuring a branch cut in the phase of the complex non-equal time correlation function g q (t), the device noise seen in Fig. 6 does not effect the identification of the step-wise change in the topological index at all. If the signal is good enough to contain the holonomy of g q (t), perfect recovery of the DQPT is guaranteed (or not at all if the quality threshold is not met).

III. ENTANGLEMENT TOMOGRAPHY
Entanglement structure and state fidelity of nonequilibrium states can be obtained following Refs. [160,[164][165][166]169], based on random measurement. Our circuit-based approach is summarized in Fig. 1(c). One can compute Rényi entropies and fidelities, and reconstruct the reduced density matrix using a generalization of the Bisognano-Wichmann (BW) theorem [180,181]. The second order bipartite Rényi entropy is where the reduced density matrix of the first N/2 sites For the quantum experiments, one may define a fidelity (of the whole system), to assess the performance of the implementation. Here, being the result extracted from experiment and ρ exact being the exact density matrix, classically computed using exact diagonalization.
To compute Eq. (22) and Eq. (23), random singlequbit basis rotations u i on each qubit i are performed using a shallow circuit, [208], followed by measurement of the probabilities P U (s) of a bit sequence s in the resulting basis. Both bipartite Rényi entropy, Eq. (22), and fidelity, Eq. (23), can be obtained using the relation [166] Tr(ρ 1 ρ 2 ) = ⟨2 N s1,s2 where D(s 1 , s 2 ) is the Hamming distance between s 1 and s 2 , and ⟨. .
. . is the average over n CUE random circuits.
Entropy and fidelity measurements are determined in position space, so the corresponding circuit involves transforming the time-evolved state from momentum to position space, as shown in Fig. 1(c) and detailed in Appendix B. The top panel of Fig. 8 shows the Rényi entropy S Note that the degeneracy between the two middle states in the spectrum is broken due to imperfections of the hardware implementation. The bottom panel shows the Bhattacharyya distance ∆B(t) between exact BW results (crosses), and simulator (blue circles) versus IonQ data (red squares).
measured in system A versus B. Both should be equal in the limit n CUE → ∞, but are not because of finite n CUE . This is an accurate estimate for the simulator data, but note that neither this difference nor the statistical shot-noise uncertainty account fully for the deviations on IonQ, which are mostly systematic.
The middle panel of Fig. 8 shows the fidelity F, extracted from the same data. Here, n shots = 1000 and n CUE = 25 result in near-perfect fidelity with the simulator. The device returns a fidelity of ∼ 50% (as the transpiled circuit consists of 26 CNOT and 84 singlequbit operations, see Table II, and hence decoherence is expected). The bottom panel of Fig. 8 shows the Rényi entropy S A+B (t) of the whole system, expected to vanish if the system is an ideal/pure state. The IonQ result deviates from zero, and shows an entropy equivalent to 1 to 2 bits of entanglement (with the environment). Statistical uncertainties are shown in panels (b) and (c) respectively, following the procedure outlined in Appendix E.
The reduced density matrix ρ A (t) of system A is reconstructed using the entanglement hamiltonian tomography (EHT) protocol of Ref. [169], based on the BW theorem [180,181] generalized to non-equilibrium states.
To do so, ρ A (t) is parametrized as [209] in terms of an entanglement hamiltonian (EH), Inspired by the BW theorem, H j are the local operators in Eq. (13) ('energy densities') and T i are commutators of the latter. In practice, at e = 0, these are with sites n, n + 1, n + 2, etc., in system A [210]. Extracting ρ A (t) from P U (s) at every t requires minimizing a functional [169] to obtain optimal parameters {β i , µ i }. All the parameters are kept independent in the fit, not assuming a specific functional form as done in Ref. [169]. The optimization is done classically and will become costly for much larger (sub-)systems. Nonetheless, scalable, efficient shadow and entanglement tomography schemes are currently under investigation, see e.g., Refs. [170,211,212]. The result of the optimization is shown in Fig. 9 for N = 4 (N A = N/2 = 2), |m|a = 0.9, and e = 0, displaying the reconstructed Schmidt spectrum p λ (t) of between the exact and BW simulator results, and between BW simulator and the BW IonQ results. The hierarchy and time dependence of p λ (t) is qualitatively reproduced, but the largest eigenvalues are systematically too small, consistent with ρ A (t) appearing more mixed.
In the presence of noise, time evolution is described by an effective but unknown Hamiltonian that is not reflected in the BW ansatz, leading to a systematic error in the reconstruction of density matrices. Errorbars are from the shot-noise uncertainty. The 'exact BW' results (black crosses) are what would be measured in the n CUE → ∞ and infinite shot limits; more discussion can be found in Appendix E. In Appendix F, we discuss device errors and use the entanglement protocol to investigate sources of decoherence. Using simple ansätze for local depolarization, bit-flip, and phase error channels, it is found that such a simple error model does not account fully for the deviations observed, and the noise mechanism in the hardware needs a more complex modeling to accurately describe the features observed in experiment [213]. In Appendix G, an EHT computation for N = 8 is presented, but the observed fidelity is essentially zero given the large number of gates used. Therefore, the entanglement structure cannot be reliably reconstructed in this case. Finally, connecting the two main aspects of this manuscript, non-equal time observables and entanglement tomography, it is shown in Appendix G how the ancilla-based interferometry scheme in Sec. II A can be avoided by reconstructing |L(t)| 2 from random measurements.

IV. CONCLUSIONS
With the advancement of quantum-computing hardware and growing access to such systems, a primary goal among domain scientists is to identify those applications of quantum computing that can still benefit from nearterm noisy and small-to intermediate-scale hardware, but are beyond capabilities of the most advanced classical computers. In the domain of lattice gauge theory, i.e., the numerical suite for simulating fundamental forces in nature such as quantum chromodynamics, computations that suffer from severe sign and signal-tonoise degradation problems appear to not be classically tractable, but they may be advanced using quantum simulators and quantum computers. Nonetheless, the size of required simulations is yet far beyond capabilities of near-term hardware. As such, it is critical to establish near-term goals of the quantum-simulation program for QCD. One such area is computing qualitative features of simpler gauge theories with similar features as QCD to discover hints of unexplored phenomena, e.g., exhibited in out of equilibrium states, that could impact our interpretation of experimental and observational data in the field of hot and dense QCD [7]. Arguably, topological phenomena with robust signatures in presence of noise and other systematics, as well as static and dynamical (quantum) phases and phase transitions, and the mechanism for equilibration and thermalization are among those near-term goals. We emphasize that, because of a typical non-logarithmic volume-law entanglement growth in these cases, most non-equilibrium phenomena are out of reach for tensor-network techniques [214][215][216][217], even for simple 1+1-dimensional models. In this manuscript, within a simple prototype model of QCD, we target quantities with such characteristics, and concretely lay out the path to obtain non-equilibrium and topological features using quantum hardware. The first results on these features are obtained by performing the simulation on commercial trapped-ion quantum hardware accessible through Cloud services. This allowed for testing the robustness of the expected features to realistic hardware noise.
The first part of this manuscript focused on non-equal time observables and correlation functions to unravel a far-from-equilibrium dynamical quantum phase transition in the lattice Schwinger model with topological features. An optimal time-evolution scheme was employed, based on a combined Fourier and Bogoliubov basis transformations, with zero Trotter error in the non-interacting limit. The scheme also allowed preparing (ground and excited) eigenstates of free-fermion Hamiltonian exactly. The algorithm can be easily generalized to a wide class of fermion models in any dimension [218], e.g., in condensed matter and atomic, molecular and optical systems. Further, a Ramsey interferometry scheme was used to quantum compute real-time observables such as Loschmidt echos and non-equal time correlation functions, at both zero and strong couplings. The underlying topological origin of the DQPT was revealed by the noise-robust holonomy of the obtained NECFs. The second part of this manuscript utilizes a random-measurement-based entanglement-tomography protocol [160][161][162][163][164][165][166][167][168][169][170][171][172][173][174][175][176][177][178][179] to extract, for the first time for a lattice gauge theory, Rényi entropies and fidelities, and, based on a generalization of the Bisognano-Wichmann theorem, reduced density matrices and entanglement Hamiltonians. These techniques are applicable to other digitized models and in higher dimensions.
The quality of results obtained using IonQ's device is promising given the indirect, and thus not fullyoptimized, access to the device. The all-to-all connectivity of the trapped-ion systems were beneficial to the scheme of this work, especially for non-equal time interferometry where a single ancillary qubit is used to control evolution on all qubits. The classical entanglement tomography studied in this work allowed the quantification, based on a few test cases, of the performance of the machine for a given gate count and entanglement of the state. Because of the indirect access, hence limited knowledge of the calibration time frame and the exact operations compiled at the hardware level, as well as trap, laser, ions, and environment characteristics, a detailed modeling of noise was not possible, and a simple error model based on single-qubit depolarization, bit-flip, and phase errors appear not to be sufficient to explain deviations from an ideal simulation. The observed deviations are likely mainly from imprecise 1-and 2-qubit gates, noting that especially the latter involves long pulse operations and are thus susceptible to device parameter drifts, changes in the position and heating of the ions, and other decoherence effects. Furthermore, the employed entanglement protocol is very sensitive to the precision of 1qubit rotations [166], and no optimization was applied in our tomography protocol against this effect. More direct access to machine parameters and control may allow for the mitigation of these effects significantly. It may be valuable to repeat and extend the presented computa-tions to other available quantum hardware in future work to cross-examine the simulations against various sources of error.
This work opens up several avenues for explorations on theoretical, algorithmic, and hardware implementation fronts. It will be interesting to extend this investigation to models with non-Abelian symmetries and to higher-dimensional gauge theories, in search for unique non-Abelian or 2D/3D features in non-equilibrium properties of matter. For realistic QCD scenarios, finite density and finite temperature must enter the considerations, and the present algorithms for non-equal time correlation functions and topological charges can be combined with efficient thermal-state preparation algorithms towards that goal, see e.g., Refs. [219,220]. Furthermore, a barely explored direction is the entanglement structure of lattice gauge theories, with important applications in the thermalization of Abelian and non-Abelian gauge theories. The reason is that entanglement structure is known to be a powerful diagnostic tool [128], and may allow characterizing quantum phases [209]. Finally, to probe entanglement structures of gauge theories on quantum computers, more efficient algorithms for entanglement tomography can be employed to take advantage of global and local symmetries of the systems [221].  In this Appendix, analytic results are obtained in the non-interacting theory, i.e., Schwinger model in the e = 0 limit, for quantities such as the Loschmidt echo, nonequal time correlation functions, Rényi entropy, and the entanglement spectrum. This is achieved by analytically diagonalizing the free fermion Hamiltonian in momentum-space computational basis. For e/|m| > 0, results are numerically determined using exact diagonalization (QuSpin [222]), as are shown in various plots throughout the main text.
For the time-evolution algorithm presented in the main text, the initial state is prepared in the (a q , b † −q ) Teigenbasis of H 0 with mass +m. This is followed by transforming into the eigenbasis of Ψ q = (ψ even q , ψ odd q ) T (the inverse transformation relative to Eq. (A4)), then transforming back into the eigenbasis of H 0 (−m). Both transformations combined lead to the 'quench' dynamics studied in this work, see Appendix A 2 for details. Hereafter, the transformation from the (a q , b † −q ) T -eigenbasis of H 0 (at either ±m) to the (ψ even q , ψ odd q ) T -eigenbasis will be called 'forward' transformation, while its inverse, from the (ψ even q , ψ odd q ) T to the (a q , b † −q ) T -eigenbasis will be called 'backward' transformation.
Note that while this work considers only two specific values of the θ angle, namely θ = 0 corresponding to mass m and θ = π corresponding to mass −m, the algorithms of this work are applicable to arbitrary values of θ. According to Eq. (3), the only change compared to the θ = 0 case is the introduction of the θ-dependent mass m cos θ, which leads to different values ω q and β, without changing the Hamiltonian diagonalization routines and subsequent analyses.

Analytics: Loschmidt echo and non-equal time correlation functions
The quench studied in the main text is realized by preparing the system in the ground state of H 0 (+m) with positive mass +m (Eq. (15)), followed by time evolution using H 0 (−m) with negative mass −m. Because the dispersion ω q is independent of the sign of the mass, H 0 (+m) and H 0 (−m) are formally identical in their respective momentum-space eigenbasis. The quench is therefore performed by explicit basis change from the H 0 (+m) to H 0 (−m) eigenbases. This transformation is q B † q (−m)B q (+m) ≡ Q, with B q (m) defined in the previous section. Time-dependent observables, after this quantum quench, can be computed analytically, which is discussed below [149].

The entanglement Hamiltonian then follows
Here, an immaterial (time-dependent) additive number has been dropped. Because H A (t) is dimensionless, this constant has no further physical meaning; conventionally the level distribution of H A (t) is scaled to lie between a standard interval, e.g., [0,1]. To obtain the entanglement spectrum, one diagonalizes H A simply by diagonalizing G A , yielding where G n are the N A eigenvalues of the (single-particle) correlation function G A (t). and ξ † n (t)ξ n (t) is the timedependent number operator in the eigenbasis of H A (Schmidt basis). Because it is in the form of a singleparticle problem, G A (t) and thus Eq. (A19) can be numerically diagonalized easily for an arbitrary large system. Computing the Schmidt spectrum p λ (t), where λ ∈ [0, 2 N A −1 ] labels the time-dependent Schmidt vectors |λ(t)⟩, or the corresponding Rényi or von Neumann entropies, from Eq. (A20) is a simple counting exercise.

Appendix B: Quantum Algorithms
The algorithmic details of the time-evolution scheme, Ramsey interferometry circuits, as well as the entanglement tomography scheme are discussed in this Appendix. All circuits are publicly available at https://gitlab. com/Niklas-Mueller1988/dqpt_2210.03089.

Time evolution
Time evolution is implemented in two bases: in the eigenmode space of H 0 in Eq. (15), and in position space, in which the interaction part, Eq. (11), can be implemented more straightforwardly. Digitization of the timeevolution operator follows a Trotter scheme, i.e., Eq. (20), which involves Bogoliubov and fermionic Fourier transforms.
The forward Bogoliubov transformation decomposes into N/2 2-qubit unitaries, B ≡ ⊗ N/4−1 q=−N/4 B q . The matrix representation for the forward transformation, in the two-qubits basis constituting momentum mode q, is where the rows correspond to the Fock basis states The action of B on these states obtains the two-qubit states ψ even q † ψ odd q † |0 even q ⟩ |0 odd q ⟩, ψ even q † |0 even q ⟩ |0 odd q ⟩, ψ odd q † |0 even q ⟩ |0 odd q ⟩, and |0 even q ⟩ |0 odd q ⟩, respectively. These states then correspond to the qubit computational basis |00⟩, |01⟩, |10⟩, |11⟩, respectively. With this convention, one arrives at the following circuit where R z/y γ ≡ e −iγσ z/y /2 with α and β defined in Appendix A 1.
The (staggered) fermionic Fourier transform (Eq. (16)), from momentum space to position space, F , is given by for the N = 8 example, where F † N/2 are separate Fourier transforms, from momentum to position space for the N/2 even and odd sites. Pauli-Z gates account for the fact that the Fourier transform is symmetric, i.e., q ∈ [− N 4 , N 4 − 1]. In the following, we will discuss the F N/2 components implementing the asymmetric Fourier transform, i.e,. q ∈ [0, N 2 − 1], keeping in mind that the aforementioned Z gates must be incorporated subsequently. fSWAP gates permute fermionic modes, accounting for their parity, so that each F N/2 has either even or odd sites as input. Examples are, in the case of N = 4 lattice sites, where rows and columns denote basis states in the following order: Here, ψ 0 and ψ 1 stand for any two fermionic modes the sub-circuit acts on, the notation indicating that qubit 0 comes before qubit 1 when circuits are read from top to bottom. fSWAP 2 , up to an overall constant phase, has the circuit representation where R x γ ≡ e −iγσ x /2 . The fermionic Fourier transform F N/2 , from position to momentum space (the inverse is used in the main text) for a staggered component, even or odd, is where n ∈ [0, N/2 − 1] labels N/2 even or odd sites and q ∈ [−N/4, N/4 − 1]. To simplify the discussion, we will from now on assume n ∈ [0, N/2 − 1] as well as q ∈ [0, N/2 − 1]. The symmetric case, i.e., q ∈ [−N/4, N/4 − 1], is realized through additional Z-gates as discussed above.
The transformation is defined recursively and consists of 2-qubit unitaries. For the simplest case N/2 = 2, the transformation is (B4) Using explicit notation, the transformation acts as follows on (two-mode) states, from (staggered) position, n to momentum, q space, where |0 q=0 ⟩|0 q=1 ⟩ and |0 n=0 ⟩|0 n=1 ⟩ are zero-fermion states in momentum and position space, respectively. Equations (B4-B5) are realized by where the row ordering is as in Eq. (B5) (from top to bottom). For N = 4, (N/2 = 2), Eq. (B6) is the full (even-or odd-site) position-to-momentum Fourier transform. Note that in our staggered formulation with N = 4 sites, this matrix is surrounded by swap gates ensuring that even (odd) sites are separately transformed, and Z gates ensuring q ∈ [−N/4, N/4 − 1]. For N = 8 (N/2 = 4), the Fourier transform of four (even or odd) modes is . A recursive pattern is evident, which can be generalized to arbitrary even N > 8. Concretely, for N = 8 and considering that W 2 4 = −1 and W 3 4 = −W 1 4 , the position-to-momentum Fourier transform F 4 is, with circuit representation where P (γ) ≡ diag(1, e iγ ) is a phase gate. This concludes the construction of the fermion Fourier transform, from position to momentum space. As part of the time-evolution scheme, first the ground state of H 0 (+m), |Ω⟩ = q |0 a q ⟩ |0 b q ⟩, is prepared in its eigenbasis. Evolving with H(−m), a change of basis between the eigenbasis of H 0 with ±m is performed, (B9) with B q andB q given by Eq. (B1). The transformation Q q has the following circuit representation, where γ 1 ≡ 2β − π 2 , γ 2 ≡ −α − π 2 , and α and β are evaluated at |m|.
For e = 0, the time-evolution circuit is (c.f. Fig. 1(a)) . This is realized, apart from a global phase, as follows where R z are single-qubit Z rotations, and [·] a/b denotes that two modes (2 qubits), a and b, are involved for every q.
In the interacting case, e/|m| > 0, following the basis change Q, the following Trotter scheme, corresponding to Eq. (20), is employed The time-evolution operator for e −iH I δt includes H I ≡ ae 2 N −1 n,m=0 ν(d nm )Q n Q m ,, where ν(d nm ) are defined in Eq. (12) and Q n ≡ ψ † n ψ n − [(−1) n − 1]/2. In position space the interaction part is diagonal and consists of 1-qubit R z and 2-qubit R zz rotations. Its implementation requires careful bookkeeping but is otherwise straightforward. Figure 1(b) summarizes the Ramsey interferometry scheme utilized in Sec. II A to compute Loschmidt echos and non-equal time correlation functions. To quantum compute Loschmidt echos, one can use

Ramsey interferometry
where measuring ⟨σ x ⟩ (⟨σ y ⟩) returns the real (imaginary) part of L(t). The controlled time-evolution operator e −iHt is either the non-interacting one, e = 0, which involves no Trotterization ('fast-forwarding'), or uses the Trotter scheme at e/|m| > 0. We emphasize that, Q, as well as (for e/|m| > 0) Bogoliubov and Fourier transforms are not controlled by the ancilla. The circuit involves controlled R z rotations for e = 0 with simple circuit implementation, and controlled R z and R zz rotations at e/|m| > 0, also with simple implementation but significantly larger circuit depth. Again, 'symm-based postselection' indicates measurements as part of an errormitigation scheme, whereby 'unphysical', i.e., particlenumber-violating measurements, are discarded. Quantum computation of the NECFs, g q (t), makes use of the decomposition where |Ω⟩ ≡ q |0 a q ⟩|0 b q ⟩. A circuit to compute both terms separately, acting on the same input states as before, is Because a q , a † q , b −q , b † −q are non-unitary, they cannot be implemented without the use of additional ancillas. This is beyond the capabilities of the device and a different strategy is employed. First, via a Jordan-Wigner transformation, consistent with the fermionic mode ordering, fermion Fock operators are mapped onto spins. Note that q ′ ; a/b indicates that the Jordan-Wigner string includes σ z q Pauli operators corresponding to both a and b modes, see the discussion of mode ordering in Sec. II of the main text. The following decomposition allows to write every term of Eq. (B11) as a sum of four unitary operators, which are separately computed. For example, The second term in Eq. (B11) is obtained in the same way. Note that the Pauli operations (and not timeevolution operators) will need to be controlled in the full Ramsey interferometry circuit shown above. Putting everything together, obtaining g q (t) involves 16 circuits, including real and imaginary parts. Because of the significant gate complexity in the interacting case, only the results for e = 0 are obtained in the main text. In this case, so that time evolution is performed mode-by-mode, each requiring two qubits and one ancilla, with a † q = σ + a,q , b † −q = −σ z a,q σ + b,−q , etc. A symmetry-based errormitigation scheme is employed, identical to that for L(t), see Appendix D for details.

Entanglement tomography
The entanglement tomography scheme is adapted from Refs. [166,169] and is summarized in Fig. 1(c). It involves random measurements via basis transformations consisting of single-qubit rotations U = ⊗ i u i , where u i is drawn from a circular unitary ensemble (CUE). In practice u i is given by , and the matrices u i are drawn randomly with the parametrization  The matrices u i are generated with nitary_grop.rvs(2) of the Python module unitary group from the package scipy.stats, based on the algorithm in Ref. [208]. The angles γ 1 , γ 2 and γ 3 in the circuit are determined by [166] The probabilities of measuring the bit string s, P U (s), in the random basis U are input for the (classical) postprocessing procedure to obtain Rényi entropies, fidelities, and the entanglement Hamiltonian via the BW ansatz, see Eq. (24) and (28) and the discussions in the main text.

Circuit construction and data generation
All circuits are implemented in Qiskit [225] using standard gates and were submitted through use of the Qiskit IonQ Provider, Qiskit ionq [226]. Table I shows the 1-and 2-qubit gate count of the subcircuit components employed in this manuscript. The number of gates refers to fundamental gates in Qiskit, obtained by iteratively decomposing the circuits using ecompose() }. 2-qubit gates count CNOT gates an 1-qubit gates are single-qubit rotations.
In contrast, Table II displays the total 1-and 2-qubit count of all quantum computations. 'Raw' refers to the circuits implemented as in Qiskit with gate count as in Table I, and 'transpiled' ('transp.') stands for the gate count of circuits that were transpiled specifically for the utilized hardware, using ranspile() }, and sen to the device through the Cloud. The gate count reflects this transpilations, with 2-qubit gates standing for CNOTs and 1-qubit gates standing for single-qubit rotations. We note that native-gate implementation, qubit-ion assignment, etc. is performed after one sends jobs to the Cloud, thus is beyond the user's control.  [100,101]; see also Refs. [227][228][229][230] for ground-laying works and pedagogical introduction of trapped-ion digital quantum computing.
The trapped-ion quantum computer consists of 13 171 Yb + ions which are aligned to form a linear crystal, suspended in a chip trap with a radial pseudopotential frequency of ≈ 3.1MHz, and spacing between the ions of about 4µm. In this setup, the 11 ions in the middle, which are more uniformly spaced, were utilized as qubits.
In the experiments of this work for N = 4 (N = 8) sites, 5 (9) of these qubits were used. The gates are implemented by using counter-propagating Raman laser beams capable of illuminating individual ions and utilizing the ionion coupling mediated by the collective motional modes. The native gate set that is physically executed on IonQ hardware is as follows [231]. The single-qubit gates are where GPI(ϕ) can be considered as a bit-flip rotation with an embedded phase and GPI2(ϕ) can be considered as an R x (π/2) gate with an embedded phase. The rotations around the Z axis of the Bloch sphere can also be performed, but these require no quantum operation and only involve a classically advancing or retarding the phase of the following operation in the circuit. The 2-qubit native gate is the Mølmer-Sørenson (MS) gate which can be expressed as   17) and Appendices A and B. Measuring σ x and σ y on the ancilla qubit yields real and imaginary part of L(t), respectively. Additionally measuring the N qubits representing the system, which are in a superposition of timeevolved and non-time-evolved states, allows us to perform a simple error-mitigation scheme: For the quench gate Q = N/4−1 q=−N/4 Q q , every Q q only mixes 0-and 2-particle states, and furthermore, H 0 preserves particle number. Therefore, the system's N qubits encoding |Ψ⟩ are in a product state, with the 2-qubits per q being in either in a 0-or 2-particle state. Measurements that yield 1-particle states indicate a bit-flip error of the device. These can simply be discarded and the remaining measurements re- weighted. The results of this procedure for L(t) at e = 0 are summarized in Fig. 10-11, demonstrating significant improvement. We also show the statistical error assuming a binomial distribution, which slightly increases with postselection as the overall statistics is reduced. For e/|m| > 0 and one Trotter step, where F is the fermionic Fourier transform and H I is the interaction part of the Hamiltonian in positions space.
The system state, |Ψ⟩, is in a superposition in position space, where both the time-evolved and non-evolved components have total particle number N/2 ('half-filling'). Again, particle-number measurement allows identifying bit-flip errors which can be discarded. Results are shown in Fig. 12, yielding no improvement. The same error-mitigation scheme is applied in Sec. II B, when computing the NECFs g q (t) (Eq. (7)) and the topological charge ν(t) (Eq. (10)). In the noninteracting limit e = 0, where g q (t) is computed modeby-mode with shallow (3-qubit) circuits, the effect of this symmetry-based error-mitigation scheme is minimal (only ≈ 5 − 15% of events contain number-violating er- rors, see Fig. 6). We note that this error-mitigation procedure was previously explored in trapped-ion computation of the Schwinger model with up to six qubits [62]. A technique based on Ref. [232] was also tried in the same work, which inserts random rotations between Trotter steps of evolution to average away the error. Nonetheless, no significant improvement was observed, pointing to the time de-correlated incoherent nature of noise in the current trapped-ion hardware. A number of different (hardwaresuitable) error-mitigation schemes are available to potentially improve our computation, see e.g., Refs. [233,234] for a recent overview of various techniques. These include randomized compiling and dynamical decoupling [235][236][237]. Randomized compiling can potentially improve our results, albeit it significantly increases simulation cost, which we therefore have not explored. Dynamical decoupling may be less profitable given the architecture and the dominant error source that we identify below to be inaccurate gate operations. The effect of gate errors could, in fact, potentially be amplified by this errormitigation scheme. As a result, the effective and computationally economic symmetry-based post-selection is the only error-mitigation scheme tried in this work.

Appendix E: Statistical Uncertainties
This section presents details of our statistical measurement-error estimation.
Only one data set per observable has been obtained in this work, hence the statistical uncertainty cannot be inferred directly by comparing many independent measurements. In part, this is because calibration drifts make it difficult to compare data sets at different times. We balanced getting data in quick succession with relatively stable systematic uncertainties, versus getting a large amount of data for statistical analysis.
Statistical uncertainties are estimated as follows. For observables based on single-qubit measurements, such as the Loschmidt echo and topological index in Secs. II A and II B, a binominal distribution is assumed, i.e. the uncertainty is ∆p = √ p 0 p 1 n shots , where p 0 (p 1 = 1 − p 0 ) is the probability of measuring 0 (1). Randommeasurement-based results, such as Renyi entropies and fidelities, and the EHT analysis are investigated as follows: Given the probabilities P U (s), n b = 100 bootstrap copies are created, i.e., one can sample from P U (s) to create additional 'measurements'. Statistical uncertainties are obtained for S (2) , F, and S A+B by computing Eq. (24) and Eqs. (22)(23) for every bootstrap sample individually. The uncertainty is the standard deviation in the limit n b → ∞, e.g., ∆F ≡ where the sum is over bootstrap samples. Shot-noise uncertainties are very small and difficult to discern in the lower two panels of Fig. 8. A close-up of the data is shown Fig. 13. The procedure is repeated for the EHT analysis. n b = 100 samples are used, except for N = 8 where due to the immense classical post-processing cost, only n b = 10 samples are tried.

Appendix F: Simple Error Models and Entanglement Tomography
Because of the indirect Cloud access to the hardware, it is difficult to identify the origin of deviation of experimental results from the theoretical expectations, which may be due to drifting machine parameters as well as ion movement and heating. In this appendix, via additional entanglement-tomography studies for N = 4 lattice sites, we address, in a simple error model, local decoherence. As will be shown, this simple model is not sufficient to describe the observed reduced fidelities. The three simple local decoherence error models employed replace the density matrix with modified ones (2) A+B = 0 is due to the finite number of random basis changes nCUE. While our allocation with IonQ did not allow us to take more data, this figure shows that the combined statistical error is very small.

as [169]
where p 1,2,3 ∈ [0, 1/N ], Tr i is the partial trace over qubit i, and N is the number of qubits. A more complete approach would be solving a Lindblad equation of the quantum circuit, see e.g., Ref. [154] for a study of the presented model, but this route will not be explored here.
The EHT analysis using the BW ansatz presented in section III is repeated using M[ρ A ] instead of ρ A in the fit, separately for three local decoherence error models Eqs. (F1-F3). Here, ρ A is given by the BW ansatz in Eq. (27) as before, and p 1,2,3 are additional free fit parameters. Figure 14 shows the results this study, using the same data as in the main text for |m|a = 0.9, e = 0, n CUE = 25, and n shots = 1000, generalized to include local depolarization via Eq. (F1). Other error sources are not considered. IonQ data without allowing for depolarization noise in the fit (red symbols) are compared to results using Eq. (F1) (green symbols). Blue symbols are simulator data and black curves denote exact results. The middle and bottom panels show Bhattacharyya distance ∆ B (t) (relative to the exact BW results) and the depolarization probability p 1 (t). The quality of the data is not improved by the ansatz and p 1 (t) varies signifi- cantly from data point to data point. Figures 15 and 16 show the same analysis for the bit-flip (Fig. 15) and phase ( Fig. 16) error channels. Likewise, including the bit-flip and phase error channels does not substantially improve the reconstruction of ρ A (t). Finally, Fig. 17 shows Rényi entropy, extracted from Eq. (24), versus that computed from ρ A (t) with the BW ansatz. Shown are fits without error channels (orange diamonds), local depolarization (green filled diamonds), local bit-flip errors (green empty diamonds), and phase errors (green empty circles). Blue symbols show the simulator data while black curves represent exact results.
Overall, these plots demonstrate that simple ansätze in Eqs. (F1-F3) fail to explain the deviations of the IonQ results. Instead, a major source of uncertainty appear to be imperfections of 1-and 2-qubit gates. Especially, 2-qubit entangling operations involve long pulse durations, and the ions may heat and change their positions as well machine parameters may drift during long circuit operations. Such effects require a more detailed modeling of the noise, which nonetheless requires more insight into the machine characteristics and the underlying operational framework that are generally not available to Cloud users. Finally, it is worth noting that the precision of 1-qubit rotations is important for the accuracy of the EHT scheme. Following the strategy of Ref. [166], improving robustness against miscalibration and drift by concatenating random unitaries did not lead to a significant improvement but will be investigated in greater detail in future work. Finally, state preparation and readout errors may need to be studied in this context too.

Appendix G: Further Results on Entanglement Tomography
The results of the entanglement tomography analysis for the largest system with N = 8 lattice sites are presented in this appendix. The significant gate count, as evident from Table II,  curves), simulator data (blue symbols), and IonQ data (red symbols). The middle and bottom panel show fidelity and total entanglement entropy, respectively. The horizontal red line indicates the expectation from a maximally mixed state, ρ A (t) = I/2 N A . Near complete failure to quantum compute the desired target state is apparent.
Regardless of hardware implementation, the BW ansatz for the reconstruction of the reduced density matrix still works reasonably well for the case of N = 8. As shown in Fig. 20, using a (perfect) hardware simulator with n CUE = 25, n shots = 2000, minimal deviations between the exact (black curves), and the BW ansatz (ideal: black plus symbols, simulator: blue circles) are observed for the reconstruction of the time-evolved state ρ A (t). While the agreement is excellent overall, these deviation can systematically be improved by including slightly less local terms in the BW ansatz in Eq. (27), as discussed in Ref. [169].
Finally, performing random measurements provides another avenue for computing non-equal time observ- The second-order Rényi entropy extracted from data via Eq. (24) (red symbols) versus that reconstructed from the BW ansatz via Eqs. (25)(26)(27)(28) with no decoherence effect included (orange diamonds), and including local depolarization, bit-flip, and phase error channels Eqs. (F1-F3) (green symbols). Blue symbols show simulator data, while black curves represent the exact results.
ables. Specifically, noting that |L(t)| 2 = |⟨ψ(0)|ψ(t)⟩| 2 = Tr[ρ(0)ρ(t)] , Loschmidt echo and rate function are computed from Eq. (24) using P U (s) obtained at different times, instead of an ancilla-based interferometry scheme. The result is shown in Fig. 21, reconstructed from the same data as in Figs. 8 and 9. The IonQ results (red symbols) are significantly less accurate than with the Ramsey scheme, because of a larger gate count compared to the direct computation in section II A.
B defined in Eq. (22)), for N = 8, |m|a = 0.8, e = 0, nCUE = 25, and n shots = 1000, including simulator (blue symbols) and IonQ (red symbols) results. The middle panel depicts fidelity F(t) (Eq. (23)). The bottom panel shows the Rényi entropy of the full system, S A+B (relative to the environment). Blue horizontal lines in the middle and bottom panels indicate ideal results, a horizontal red line indicates zero fidelity or maximal entropy (ρ(t) = I/2 N ), respectively. See Fig. 19 for a close-up of the lower two panels.