Real- and imaginary-time evolution with compressed quantum circuits

The current generation of noisy intermediate scale quantum computers introduces new opportunities to study quantum many-body systems. In this paper, we show that quantum circuits can provide a dramatically more efficient representation than current classical numerics of the quantum states generated under non-equilibrium quantum dynamics. For quantum circuits, we perform both real- and imaginary-time evolution using an optimization algorithm that is feasible on near-term quantum computers. We benchmark the algorithms by finding the ground state and simulating a global quench of the transverse field Ising model with a longitudinal field on a classical computer. Furthermore, we implement (classically optimized) gates on a quantum processing unit and demonstrate that our algorithm effectively captures real time evolution.

The current generation of noisy intermediate scale quantum computers introduces new opportunities to study quantum many-body systems. In this paper, we show that quantum circuits can provide a dramatically more efficient representation than current classical numerics of the quantum states generated under non-equilibrium quantum dynamics. For quantum circuits, we perform both real-and imaginary-time evolution using an optimization algorithm that is feasible on near-term quantum computers. We benchmark the algorithms by finding the ground state and simulating a global quench of the transverse field Ising model with a longitudinal field on a classical computer. Furthermore, we implement (classically optimized) gates on a quantum processing unit and demonstrate that our algorithm effectively captures real time evolution.

I. INTRODUCTION
Ground states of strongly correlated systems and their quantum dynamics far from equilibrium present important problems in understanding quantum matter. In both cases, we often rely upon numerical tools to unravel the emergent physics. Our most general tool is exact diagonalization (ED), which is limited in its scope because it requires storing an exponential number of parameters with respect to the system size [1]. Besides ED, one can efficiently find the ground states or simulate dynamics of one-dimensional local gapped Hamiltonians using matrix-product state (MPS) techniques, such as the density matrix renormalization group (DMRG) algorithm and the time evolving block decimation (TEBD) algorithm [2,3]. However, for generic systems, the rapid growth of entanglement under far-from-equilibrium dynamics severely limits the accessible time scales due to the cost of storing or sampling the state. For systems without the infamous sign problem, quantum Monte Carlo techniques represent a powerful tool [1,4]. Importantly, many physically interesting systems fall outside the scope of these modern numerical methods and new approaches are needed to tackle these.
Universal quantum computers have become an increasingly feasible setting for simulating quantum dynamics [5,6]. Current Noisy Intermediate Scale Quantum (NISQ) devices contain of order 50 qubits and give access to hundreds of quantum gate operations [7]. NISQ devices have a fundamental advantage over classical numerics -the physical resources required to store quantum states grow linearly, not exponentially, with the system size. Although the noise precludes implementing many quantum algorithms, it is believed that the simulation of quantum systems and dynamics may be one of the most powerful uses of NISQ quantum computers before scalable error correction is implemented. Indeed, there have already been several works demonstrating the use of quantum computers for this purpose that have benchmarked the currently available devices [8][9][10]. These works demonstrate that it may be possible to study classically inaccessible systems on near-term generations of quantum computers. Algorithms that offer methods to study quantum systems using NISQ devices are thus of significant interest.
Experimental advances in quantum computation technology have also raised several fundamental questions about the relationship between complexity and entanglement of physically relevant quantum states [11]. In classical algorithms, especially tensor network methods, the entanglement is a good proxy for the difficulty of representing a state. For a quantum circuit, however, these measures are relatively independent; one can have states with high entanglement but low complexity. This distinction between complexity and entanglement means there is a complexity window, between those states with high entanglement that are accessible with polynomial depth circuits and those with circuit depth exponential in system size [12], shown schematically in Fig. 1(b). While the former marks the limit of current classical numerical methods, quantum simulators and computers may allow us to study a new class of physically interesting states in this complexity window.
In this paper we study a class of quantum circuits motivated by a representation of matrix-product states. For a given amount of entanglement, this class requires exponentially fewer parameters. The structure of this paper is as follows. we demonstrate that the ansatz states are a good approximation for the states obtained during time evolution in Section II by comparing with classical numerics. In Section III, we consider a variational time evolution algorithm for general quantum circuits for both real and imaginary time evolution. In Section III C, we classically optimize the gates, then implement the compression directly on a quantum processing unit (QPU). We conclude, in Section IV, by noting several future avenues of exploration using the techniques developed in this work. The space we capture can be generically characterized as states with low complexity but high entanglement, which efficiently captures time evolution. (c) To perform time evolution, we prepare a state |Ψ M qc (t) , then apply a Trotterized time evolution to obtain |Ψ M qc (t + ∆t) . By variationally optimizing each of the gates, we find an optimal representation of the time evolved state within the sub-manifold defined by our variational ansatz.

II. COMPRESSED CIRCUITS
Although entanglement is a good proxy for the difficulty of representing a state using an MPS ansatz, the light cone determined by a time evolution under a local Hamiltonian enforces a particularly simple entanglement pattern that can in principle be captured with fewer parameters. We use an ansatz where sequential quantum circuits represent our states. These circuits consist of a set of two qubit gates {U i } that are applied sequentially as shown in Fig. 1(a). The circuit is said to be of order M when there are M "layers" of gates. The total depth of this circuit is 2(M − 1) + N − 1, which scales linearly in both the system size N and with the order M . In contrast to a more commonly studied brickwall circuit structure [13], this ansatz does not restrict the correlation length of the states we can represent, see Appendix A for additional discussion.
We note that the states defined by these quantum circuits form a sub-manifold of matrix-product states with bond dimension χ = 2 M . In the case of M = 1, the quantum circuit is exactly equivalent to an MPS of bond dimension χ = 2 (see Appendix A). However, for M > 1, these quantum circuits have exponentially fewer param-eters than a generic matrix-product state in canonical form with bond dimension 2 M . In other words, these quantum circuits describe states with high entanglement but low complexity, which-as we demonstrate belowencompass time evolved states. Note that this reduction of parameters does not necessarily translate into a sparse representation when stored as an MPS on a classical computer.
To test this class of quantum circuit ansatz, we first consider far-from-equilibrium dynamics of a global quantum quench. Crucially, such dynamics is typically accompanied by fast ballistic growth of entanglement, which puts mid-to-long time dynamics out of reach for numerics beyond small systems. Concretely, we consider dynamics under the Hamiltonian which is a quantum Ising spin chain on N sites with both transverse (g) and longitudinal (h) fields. For the special case h = 0, the model is integrable. We consider a global quantum quench protocol with polarized initial state |Ψ = | · · · ↑↑↑ · · · at time t = 0. Our goal is then to accurately approximate the state |Ψ(t) = e −iĤt |Ψ , at (real or imaginary) time t after the quantum quench.

A. Efficient Representation of Quantum States
We now demonstrate the representation power of the quantum circuit ansatz by comparing it with classical numerics using MPS. We first perform the time evolution using 4 th order Trotterized TEBD [3] for N = 31 with maximum bond dimension χ = 1024 and step size τ = 0.01 to obtain quasi-exact approximation of the state |Ψ(t) . This bond dimension ensures that our results are close to exact for all considered timescales. We then take the MPS at a selection of times, which we denote |Ψ mps (t) , and find the optimal quantum circuit of order M , which we denote |Ψ M qc (t) . The state represented by the quantum circuit is implicitly parameterized by a set of two-qubit unitaries {U i (t)}. We perform an optimization over the unitaries in our quantum circuit to find the state with maximum fidelity This is done by iteratively by updating each U i (t) using a polar decomposition [14] (see Appendix B for more details).
In Fig. 2 we show the fidelity of the quantum state obtained from the quantum circuit ansatz as well as the half-chain von Neumann entanglement entropy, S. are chosen such that the dynamics of the system are expected to be chaotic and hard to simulate due to fast scrambling [15,16]. The accuracy of the approximation decreases with time as correlations build throughout the system, but improves as the order M is increased. For a given order M , this data also shows that the circuit more accurately captures the state for weaker h, indicating an increase in the complexity of the simulation for larger h. Figure 2 also shows the growth of entanglement. We find that the ansatz easily captures the rapid ballistic growth of entanglement for small h. As we increase h, we find that the growth of entanglement slows down. This indicates that the practical complexity of the quantum states increases with h whereas the growth of entanglement decreases, thus partially closing the still exponentially large complexity window.
From this data we can compare the number of parameters required to achieve a given accuracy using our quantum circuits with those needed for an MPS. For a given order M we find the time t * up to which the fidelity is greater than F = 1 − 10 −4 , indicated by the grey dashed line in Fig. 2. In Fig. 3, we plot the number of parameters in the quantum circuits and the MPS as a function of the reachable time t * . This figure shows that the number of parameters in our quantum circuit ansatz scales linearly with the reachable time t * , in stark contrast to the exponential growth in parameters for the MPS. Note that the circuit depth of a fully Trotterized time-evolution also scales linearly with time [17,18] and has for sufficiently small time steps an error for local observables that is independent of system size as well as simulation time [19]. However, we find that the compressed circuit generically performs better while the quantitative improvement over the fully Trotterized time-evolution depends on the model parameters-this reduction of circuit depth is particularly valuable for current NISQ devices on which Trotterized time evolution is very challenging [10].
We stress that the linear scaling of the number of parameters persists across the different values of h. The additional complexity for large values of h appears as a change in the gradient of the linear scaling. These results demonstrate that the complexity of the quantum state grows linearly in time, while the MPS ansatz requires a number of parameters that grows exponentially in time due to the linear growth of entanglement. For all values of h we can see that the quantum circuit has an exponential advantage over MPS in terms of the number of parameters required. Even for short times of O(1) in the coupling J, we require fewer parameters to accurately represent the state with a quantum circuit than with an MPS.

III. VARIATIONAL TIME EVOLUTION ALGORITHM
Having confirmed the representation power of our ansatz, we now demonstrate how to implement time evolution restricted to the states defined by our ansatz. This, in turn, demonstrates that the optimization of the quantum circuit can be performed on a quantum device using hybrid quantum optimization algorithms. This potentially enables the simulation of dynamics beyond the reach of classical numerical methods, which are limited by the cost of storing the quantum state.
Our algorithm for time evolution is shown schematically in Fig. 1(c). Given the quantum circuit at time t, we apply a second-order Trotterized approximationV (∆t) = e −iĤeven∆t/2 e −iĤ odd ∆t e −iĤeven∆t/2 to the time evolution operator e −iĤ∆t . We then find the state |Ψ M qc (t + ∆t) that maximizes fidelity That is, we iteratively optimize over the set of 2site unitary gates {U i (t + ∆t)} that define the state |Ψ M qc (t + ∆t) . We perform this optimization similarly to that in the previous section, where we update each U i iteratively using a polar decomposition. After multiple sweeps, we find that the fidelity F in Eq. (3) converges. Some variant of the optimization algorithms that have already been tested on quantum devices for variational quantum eigensolvers could also be applied effectively for our algorithm [20,21], something we leave for future work.

A. Real Time Evolution
In Fig. 4 we show the local magnetization and the halfchain entanglement entropy simulated using our quantum time evolution algorithm. We consider the same quantum quench protocol as above with h = 0.1. Importantly, this case is non-integrable and has a fast linear growth of entanglement under the non-equilibrium dynamics.
Our results show that we are able to accurately capture the magnetization for times that scale linearly with the order M . Here it is important to note that the time evolution is performed entirely within the sub-manifold of circuits defined by our ansatz with a fixed order. We additionally find that we are able to capture the linear growth of entanglement using these quantum circuits and that the saturation of the entanglement depends linearly on the order M . In contrast, the corresponding MPS representation has an exponentially large bond dimension requiring O(2 M ) parameters.
We emphasize that this time evolution algorithm is different from the time-dependent variational principle (TDVP) algorithm simulating time evolution with a quantum circuit proposed in [22,23]. In those approaches, one solves the TDVP equations approximately by stochastic sampling, i.e., measurement, and performs finite time stepping by numerical integration [24]. In the present algorithm, we first perform finite time stepping by Trotterization, and then try to find the optimal states within the sub-manifold defined by our ansatz. This is much closer to the tDMRG [25,26] or TEBD [27,28] algorithms, but also has similarities with the iTDVPinspired algorithm proposed in [29]. This optimization algorithm for the time evolution is similar to one used for the multi-scale entanglement renormalization ansatz (MERA) [30], and in the context of symmetry-preserving ansätze [31].
The problem of efficiently optimizing a variational ansatz is one that is common to many current hybrid quantum-classical algorithms [32,33]. While we have classically performed the optimization by using polar decomposition, this could be replaced by any (stochastic) gradient descent based optimization scheme. This optimization of a large number of parameters within a restricted manifold is often plagued by barren plateaux [34]. Further work is therefore required to guarantee the efficiency of optimization for this wide class of algorithms, including our own.

B. Imaginary Time Evolution
We can also apply this time evolution algorithm to find ground states using imaginary-time evolution. In this section, we first explicitly show how one can embed the required non-unitary operators in unitary gates using an ancilla qubit and post-selection [35]. Second, we demonstrate that our ansatz can effectively converge to the ground state under imaginary-time evolution.
One can formally write down the exact imaginary-time evolution procedure |GS = lim τ →∞ e −Ĥτ |ψ 0 , where τ is real. This is equivalent to evolving in imaginary time (t → −iτ ) and corresponds to acting on the state with a non-unitary operator, which becomes a projector onto the ground state in the limit τ → ∞. We perform imaginary-time evolution analogously to our realtime evolution algorithm, where we sequentially compress the state back onto our ansatz as in Eq. (3) but witĥ V (∆t) = e −Ĥ∆t . Similarly to real-time evolution,V (∆t) can be approximated by a product of 2-qubit non-unitary gates using Trotterization.
To perform imaginary-time evolution, we are therefore required to implement non-unitary gates on the quantum computer. We achieve this by embedding the non-unitary gate in a unitary gate acting on one extra ancilla qubit. For a generic non-unitary operator A acting on N qubits, we define a unitary (N + 1)-unitary V A by The strategy we employ is to find a block C and a scaling factor s that ensures the first 2 N columns of V A are mutually orthonormal, which guarantees unitarity. The remaining columns can be fixed using a QR decomposition. We explicitly show the full embedding procedure in Appendix C. Note that if we were to implement a full To find the optimal performance of our ansatz, we perform a procedure similar to VQE, where we iteratively optimize the expectation value of the Hamiltonian in Eq. 1 for our ansatz. For these depths, our imaginary-time evolution algorithm successfully converges to the optimal point, depicted by the dashed lines. The χ = 4 line indicates the ground state energy of an MPS with bond dimension 4 found using DMRG. To achieve a better accuracy we successively decrease the time step ∆τ after achieving convergence for the previous step size.
Trotter step for each optimization step, as we did previously for real-time evolution, we would require a linear number of ancilla qubits resulting in an exponential cost due to post-selection. Instead, one should apply and optimise the state for each 2-qubit gate in the Trotterized time evolution separately. In this case, the total number of measurements required across all Trotterized gates in a single time step scales only linearly with system size.
To benchmark how well our ansatz can approximate the true ground state, we directly minimize the energy This procedure is similar to a variational quantum eigensolver (VQE) [32], where the parameters encoding the quantum state are iteratively adjusted to minimize the energy. We perform the procedure on a classical computer where it is intended to benchmark our imaginarytime evolution algorithm, where we consider the energy in Eq. (5) to be the best achievable by our chosen ansatz, shown as dashed lines in Fig. 5. Instead of using gradient descent methods, we iteratively replace the unitaries using polar decomposition, see appendix B.
In Fig 5, we show the results of our imaginary time evolution. These show that we can successfully converge to the optimal energy attainable with this ansatz. As expected, the results for M = 1 match those from DMRG with bond dimension χ = 2 due to the equivalence between the circuit and MPS representation. Note that while a modest MPS bond dimension χ = 4 performs better than our ansatz for M = 2, 3, we still achieve errors well below the threshold of current NISQ hardware, which validates this approach as a method for finding ground states on a quantum device.
We note that an alternative approach to imaginarytime evolution was taken in the QITE algorithm [36,37].
There it was noted that if enough information about the initial state is known, a non-unitary gate can be replaced by a unitary one without the use of ancillas. However, to get closer to the ground state requires state-dependent unitary operators with increasingly large support. In contrast, our algorithm requires a fixed set of local gates that can be repeatedly applied to reach later times, much like the TEBD algorithm for MPS [27]. While the approximation step is stochastic on a quantum computer, the overall procedure deterministically converges to the ground state. The choice of ansatz is also completely flexible. Viewing the procedure as a sequential compression in this way raises an interesting comparison with the compression of a tensor network to form a MERA and the emerging view of learning with tensor networks as a procedure of compression [38,39].

C. Simulation on QPU
While our algorithms are designed for near-term quantum computers, the noise and coherence times of currently available devices place strong limits on what can be achieved. However, we are able to demonstrate parts of the algorithm on a quantum computer by delegating more of the algorithm to the classical computer. Here we classically optimize the time evolved states |Ψ M qc (t) , then construct and measure the corresponding state on a QPU, namely the 5 qubit IBM-Q device codenamed Bogota [40,41]. This process allows us to access times on the QPU that are inaccessible using standard Trotterized evolution techniques.
Concretely, we consider the following quantum quench setup on N = 5 qubits. We initialize the system in the product state | − − + ++ , i.e. a domain wall in the x-basis, and evolve with the Hamiltonian (1) with g = 0.25, h = 0.2. For this range of parameters and initial state the dynamics is dominated by the motion of a single mobile domain wall and so can be well approximated by an order M = 1 circuit. The longitudinal field, h, leads to a linearly confining potential between domain walls, and in the case of a single domain wall corresponds to a linear background potential leading to Wannier-Stark localization [42]. In Fig. 6(a), our ED results show the characteristic periodic melting and revival of the domain wall.
In Fig. 6 we show the results of constructing and measuring our compressed quantum state on the IBM QPU compared with ED results. Here we optimize the set of gates {U i (t)} on a classical computer, which is then fed to the QPU to create the quantum state. The measurement of the magnetization in the x-basis closely matches the exact results. In particular, the spatial distribution of the magnetization (Fig. 6(a)) show the periodic spreading and reconstitution of the domain wall. Furthermore, the magnetization on the central spin, shown in Fig. 6(b) accurately and quantitatively matches the ED simulation for long times, which are not limited to the range we have considered. These timescales are currently inaccessible using a naive Trotterized evolution on this quantum device, which would require a circuit depth of O(t).

IV. DISCUSSION
In this paper, we have shown that physically relevant quantum states, namely ground states and those arising under non-equilibrium dynamics, can be efficiently represented using a sequential quantum circuit ansatz. This ansatz describes a "sparse" representation spanning a corner of the larger MPS manifold. For time evolution, the time scales that we can reach scale linearly with the number of parameters in the circuit, representing an exponential advantage over existing classical methods. This suggests that even within the class of MPS defined by a fixed bond dimension, there exist a range of physical states for which our ansatz is a more efficient representation than an MPS. To exploit the representation power, we used a time evolution algorithm for a general quantum circuit ansatz that can be implemented natively on existing quantum computers. Importantly, the quantum circuit ansatz is flexible and is not restricted to the one used in this paper. Using near term devices this may provide access to non-equilibrium dynamics beyond the reach of current classical algorithms. Finally, we have shown that this time evolution algorithm can also be applied in imaginary-time to obtain ground states on a quantum computer.
The optimization procedure that we used [14]-fidelity maximization using a polar decomposition-may have other potential applications. For instance, instead of considering the compression of states, one can consider the compression of unitaries. This technique can be applied to approximate a multi-qubit unitary by a series of 2qubit unitaries, or to compress a deep quantum circuit. Both of these are particularly important for current NISQ devices.
Our procedure is also a potential practical tool for studying quantum complexity. Quantum state complexity is an intriguing research field, but is difficult to study numerically. Previous results primarily focus on noninteracting systems [43][44][45]. By using states acquired from procedures such as TEBD and DMRG, and approximating them using a chosen ansatz and polar decomposition methods, one can concretely probe the complexity of generic classes of states (such as quantum scar states and many-body localized states) that were previously difficult to analyze. Additionally, the window between complexity and entanglement is of significant interest. In particular, Ref. [12] uses a random unitary circuit model for time evolution to demonstrate that even when the growth of entanglement saturates for a finite system, the complexity of the quantum states continues to grow linearly in time over far longer time-scales. This highlights a large window between maximally entangled states and maximally complex states. Our work shows that this window appears to shrink for non-integrable systems (see Fig 2). The techniques developed in this paper open the opportunity to directly study complexity windows in concrete systems.
The algorithms studied in this work open up several intriguing generalizations. First, one could apply the algorithm to study short time dynamics for higher dimensional systems, which are generally difficult problems for classical numerics. Applied directly on a quantum computer, this algorithm offers a tractable way to study higher dimensional systems at large system sizes and to probe physics that only manifests at higher dimensions. Moreover, the algorithms considered are agnostic to the specific ansatz used. It is an interesting question to compare how an ansatz with a different entanglement pattern performs. For instance, in Ref. [46] quantum circuits containing entangling gates acting over the full system are considered and optimized to represent time evolved states by a reinforcement learning approach, which is complementary to our approach. Additionally, quantum circuits inspired by matrix product states have shown promise for solving non-linear Schrödinger equations [47,48]. Similar analyses for various ansatz structures could shed light on the deeper relationship between entanglement and complexity. In this section, we describe an exact mapping between an MPS of bond dimension χ and a sequential quantum circuit with (n + 1)-site unitaries, where n = log 2 χ. Given such an exact equivalence, one can approximate (n+1)-site unitaries with 2-site unitaries to arbitrary precision. This results in the ansatz we consider in the main text, which corresponds to the "sparse" matrix-product states. More generally speaking, one can always rewrite isometric tensor network states as quantum circuits [49].
An MPS in right canonical form is given as, indicates each individual tensor B [k] is an isometry mapping from |α k−1 → |α k , i k . Any isometry can always be rewritten as a unitary acting on a normalized state |0 k , i.e.
where the state |0 k would have dimension dim(|0 k ) = χ k × dim(|i k )/χ k−1 . We assume dim(|0 k ) to be an integer without loss of generality since we can always enlarge the bond dimension to match this condition. One can easily verify the equivalence by substituting Eq. A4 into Eq. A2. Once the connection between the isometries B [k] and unitaries U [k] acting on a state |0 k is established, we can rewrite the right canonical MPS as a quantum circuit with a set of corresponding gates {U [k] } acting on the initial state |0 ⊗N (see Fig. 7). As expected, the dimension of the final state is the same as the initial state, because the virtual indices {α} are internally contracted.
For spin-1/2 systems, the physical dimension is d = 2 and we have a standard quantum circuit operating with qubits. If the MPS consists of tensors B [k] with bond dimension χ = dim(α k−1 ) = dim(α k ) = 2 n , where n ∈ N, then the corresponding unitaries act on (n + 1) qubits. As a result, an MPS with maximum bond dimension χ is equivalent to a quantum circuit defined by unitaries acting on maximally log 2 χ + 1 sites sequentially. These unitaries can then be further decomposed into a series of sequential 2-site unitaries where the number of required 2-site unitaries scales polylogarithmically with respect to the inverse of the desired error.
Moreover, an MPS of bond dimension χ = 2 maps exactly to our circuit ansatz of order-1. Note that this is particular to our ansatz; the commonly studied brickwall circuit structure with two layers, which has the same number of two-site gates, can be mapped to an MPS of bond dimension χ = 2 but can only represent states with finite correlation length. The above connection between matrix-product states and quantum circuits is well-known in the community and was recently applied in several works [13,29,50]. Repeating such replacement, one arrives at a circuit with pattern as in Fig. 1 (a).

Appendix B: Classical simulation algorithm for quantum circuit
In this section, we describe two algorithms. The first algorithm maximizes the fidelity between two states defined by a set of unitaries, similar to the known Evenbly-Vidal algorithm [14]. The second algorithm uses the first algorithm to perform time evolution restricted to the space defined by the ansatz under consideration.
To maximize the fidelity F = | Ψ target |Ψ M qc | 2 , we iteratively optimize the fidelity with respect to each gate U i,j , while keeping the remaining gates fixed. Note that the double indices (i, j) refer to order and site respectively, whereas in the main text we group the indices into a single index.
We first rewrite the overlap between the target state |Ψ target and the order-M circuit |Ψ M qc in the following form, where U i,j is the unitary to optimize and E is the environment matrix as shown in Fig. 9(a).
is equal to the square of the real part of the overlap. This is because any global phase offset can always be compensated by absorbing a single site rotation into the 2-site unitary. The solution to the unitary maximizing Re [ φ|U i,j |ψ ] is known; for E = XΣY † , the optimal U i,j is given by Y X † as in Fig. 9(b).
To obtain the optimal circuit, we iterate through all of the gates and update each gate with the exact solution of the local optimization problem. Given a maximal iteration number N iter , absolute convergence error a , and relative convergence error r , the algorithm is described in Alg. 1.

Algorithm 1: Maximizing overlap
Input : |Ψtarget , |Ψ M qc , Niter, a, r Output: A set of {Ui,j} maximizing Ψtarget|Ψ M qc ({Ui,j}) , idx = 0, 0 = inf; while idx < Niter and idx > a and ∆ > r do idx = idx + 1; We used standard tensor network techniques to construct the environment tensor and truncated singular values less than 10 −14 . The algorithm was made significantly less expensive by caching and updating the environments to avoid recomputing the entire environment from scratch during each new iteration. For our computations, N iter = 10 5 , a = 10 −12 , r = 10 −4 .
We now introduce our second algorithm, which performs time evolution directly on the manifold defined by our ansatz. To time evolve a state |Ψ(t) , we maximize the fidelity F = | Ψ(t + ∆t)|V (∆t)|Ψ(t) | 2 , where our unitaries parameterize |Ψ(t + ∆t) andV (∆t) is a single Trotterized time step. In this way, we can iteratively evolve forward in time from an initial state. The overall algorithm for time evolution is given as in Alg. There are two primary sources of error in our algorithm: the Trotterization error and the projection error. The Trotterization error arises from approximating the true time evolution operator by a series of 2-site gates. This can be made arbitrarily small by decreasing ∆t or by taking higher order Trotter decompositions. The projection error arises from projecting the time evolved state back onto the manifold of circuits of order M . This error is affected by the chosen ansatz and limits the time to which one can simulate within a given error threshold.
We can estimate the total error by monitoring the fidelity at the end of each optimization i F i . This total error estimate is accurate as long as the Trotterization error remains small and if F i is close to 1 at each step. As an example, in Fig. 10 we show the error estimates for the simulation performed in Fig. 4. We see that the time when the error crosses the threshold matches the time when σ z starts to deviate.

Appendix C: Non-unitary gates
In this section we describe a procedure to embed an arbitrary non-unitary N -qubit operator A in an (N + 1)qubit unitary gate. We consider the first of these (N +1)qubits to be an ancilla qubit that we initialize in the |0 state and project into the state |0 by post-selection. Our claim is that there exists a unitary of the form where s −2 is the maximum eigenvalue of A † A (or equivalently AA † . We note that the matrices A † A and AA † have real and non-negative spectra. This follows from a singular value decomposition, i.e. A = U ΣV † with U, V unitary and Σ non-negative real and diagonal, and so Our goal is to show that for any A ( = 0, although this case can also be included) we can find the 2 N × 2 N matrices, B, C and D such that U A is unitary. Our approach is the following. We first note that U A being unitary is equivalent to the statement that the columns of U A form an orthonormal basis of C 2 N +1 . We will then use this to find a block C and a scaling factor s consistent with this, i.e. such that the first 2 N columns of U A are orthonormal. Given A, C, and s we can then use a QR-decomposition to easily find B and D, as explained below.
Let us denote the columns of A and C by a j and c j respectively, e.g., [C] ij = [c j ] i . For U A to be unitary C and s must satisfy In terms of the column vectors, these can be written as For i = j this is a statement that the first 2 N columns of U A are normalized, and for i = j it is the statement that these columns are mutually orthogonal. Next we note that if C satisfies Eq. (C2), then we have the singular value decomposition C = UΣV † , where U and V are the same unitaries as in the SVD of A = U ΣV † . This implies that Since Σ 2 must be non-negative, we only have a solution to Eq. (C2) if s −2 is greater than the largest eigenvalue of A † A (all of which are non-negative), and so we set s −2 equal to the largest eigenvalue. We therefore have that Σ 2 = 1 − s 2 Σ 2 , with our choice of s ensuring thatΣ is real and non-negative. Finally, given A and C, we can find the blocks B and D using QR-decomposition. Namely, let us construct the matrixŨ whereB andD are random matrices, then by QR-decompositionŨ where U A is the unitary in Eq. (C2) and R is an upper triangular matrix. Since the first 2 N columns of U A are orthonormal they will be untouched by the QRdecomposition algorithm. gates are decomposed into a series of finitely-many gates selected from some universal gate set. A small perturbation of a two-site gate may lead to a large perturbation in the decomposition. These differences translate into large fluctuations in the measured observables due to the imperfections in the QPU.
To compensate for this problem, we average over the gauge freedom in a quantum circuit. Given the two-site gates {U i } describing the quantum states, there are gauge degrees of freedom to insert identities described by random unitaries and their complex conjugates. For example, if U i U i+1 act consecutively on the same qubit, we can insert the random single-site unitary V and its complex conjugate as and obtain the two-site gates W i+1 W i describing the same operation. To average over the gauge degrees of freedom, we average measurement outcomes corresponding to circuits differing by the insertion of random unitaries and their conjugates. This procedure mitigates the previously mentioned error to a certain extent.