Noise-robust exploration of many-body quantum states on near-term quantum devices

We describe a resource-ef ﬁ cient approach to studying many-body quantum states on noisy, intermediate-scale quantum devices. We employ a sequential generation model that allows us to bound the range of correlations in the resulting many-body quantum states. From this, we characterize situations where the estimation of local observables does not require the preparation of the entire state. Instead smaller patches of the state can be generated from which the observables can be estimated. This can potentially reduce circuit size and number of qubits for the computation of physical properties of the states. Moreover, we show that the effect of noise decreases along the computation. Our results apply to a broad class of widely studied tensor network states and can be directly applied to near-term implementations of variational quantum algorithms.


INTRODUCTION
Quantum computers offer computational power fundamentally different from classical computers. A universal quantum computer may solve classically intractable problems within areas ranging from many-body physics to quantum chemistry 1 . There has been impressive experimental progress in developing quantum computers on different architectures [2][3][4] . Although achieving faulttolerant computation remains a challenge, noisy intermediatescale quantum (NISQ) devices are expected to be available in the near future 5 . These are devices containing a few hundred qubits with small error rates but without error-correction. An outstanding question is what kind of computations such devices may facilitate.
Algorithms designed for NISQ devices should run on a moderate number of qubits and be resilient to noise. The specific hardware may also pose further restrictions regarding the connectivity of the device, as not all qubits can interact directly with each other 2,4 . Promising frameworks that fulfill these conditions are the quantum approximate optimization algorithm 6 and the quantum variational eigensolver (VQE) 7,8 . In these frameworks, the task of the quantum computer is roughly speaking to compute the expectation value of local Hamiltonians on some many-body quantum state. Recent work has characterized a number of conditions for which this can be done in a noiserobust way [9][10][11][12][13] . Due to the limited resources of NISQ devices, it is also important to run such algorithms as efficiently as possible in terms of circuit size and number of qubits.
We address this outstanding problem by developing a general framework for computing physical properties of quantum manybody states efficiently on NISQ devices. In particular, we upper bound the circuit size and number of qubits necessary to estimate the expectation values of local observables. Importantly, these bounds can significantly decrease the resource requirements compared to previous works for a number of circuit topologies and sizes. Specifically, we are able to show that the energy of a many-body quantum state can be estimated with a constant-sized quantum circuit if the correlation functions exhibit an exponential decay. This is the case for non-trivial states such as ground states of gapped Hamiltonians, surface codes and quantum states described by a multiscale entanglement renormalization ansatz (MERA) or the larger class of deep MERA (DMERA) 12,13 . The latter is believed to capture Chern insulators.
Our framework is akin to sequentially generated [14][15][16] or finitely correlated states 17 . This enables us to control the size of the past causal cone 18-20 of local observables. Combined with the notion of mixing rate of local observables under the circuit 12,13 we determine after how many layers of the circuit, expectation values stabilize. To estimate these expectation values, it suffices to implement the potentially small subset of the circuit under which they stabilize instead of producing the entire many-body state or its past causal cone. Consequently, the necessary number of qubits and quantum gates can be reduced significantly from scaling with the size of the many-body state to even a constant number.

Basic setup
We consider three basic operations, which are iterated T times to generate a many-body quantum state. The first operation adds qubits to the existing system. The second operation lets them interact with each other and a bath via a constant depth circuit. In the third operation, some of the existing qubits may be discarded. By introducing a separate bath, we allow for situations where a fixed sized quantum processor (the bath) iteratively prepares a quantum state on the system qubits. This allows e.g. the construction of arbitrary matrix product states (MPS) (see Example 1 below).
More specifically, we will start with a system S 0 consisting of n 0 qubits initialized in some fixed state ρ 0 and a bath system B consisting of s B qubits initialized in some fixed state ρ B . At each iteration t, we introduce new subsystems S t with n t qubits and ancillary states A t with a t qubits, all initialized in some fixed quantum state. These new subsystems then interact with the existing ones and finally, the ancillary system is discarded, which concludes the iteration. The procedure is iterated for a total of T iterations to produce the entire state.

Interaction scheme
The structure of the final quantum state is determined by the allowed interactions between qubits during the iterative preparation. In order to get a handle on how the correlations of the final state evolve during the preparation, we will fix the allowed interactions according to a given interaction scheme. Our construction of such an interaction scheme is inspired by socalled graph states 21 . These can be visualized by letting the vertices of a graph denote qubits and edges denoting correlations between them. In a similar way, we define an interaction scheme as a sequence of T graphs {G t = (V t , E t )} where a vertex (V t ) denotes a qubit and an edge (E t ) implies that it is possible to implement unitary gates between these two qubits in that time step. We further restrict this unitary gate such that at most D rounds of twoqubit operations are applied with the condition that on each round at most one unitary acts on each qubit. Figure 1 showcases a simple example of one iteration of the procedure illustrating all aspects of the framework. While such a simple interaction scheme serves the purpose of illustration, our framework encompasses arbitrary interaction graphs in each iteration step. This allows it to capture a number of widely studied tensor networks states, as illustrated in the two examples below.
Definition 1 (MPS and higher dimensional versions) MPS have been studied in a sequential interaction picture 22,23 and adapt naturally to our framework. The initial system S 0 consists of one qubit and the dimension of the bath gives the bond dimension.
At each iteration t, we add a system consisting of one qubit S t to the system. The graph G t only has edges between the qubits in the bath and the newest added qubit S t . Note that we are considering a proper subset of MPS since we restrict to unitaries implementable with D two-qubit gates for each edge. Nevertheless, the bond dimension of the resulting MPS scales exponentially with the number of qubits in the bath. It is straightforward to generalize such sequentially generated states by considering the case in which a bath interacts with a subsystem of dimension d rather than a single qubit at each iteration 12,15 .
Definition 2 (Deep multiscale entanglement renormal-ization ansatz) The DMERA, introduced in ref. 13 , is a variation of the MERA 19 tailored for NISQ devices. In our framework, the initial system S 0 consists of one qubit and there is no bath. We then define the graphs G t recursively: at each iteration we add one Fig. 1 Generation procedure. Example of one iteration of the generation procedure broken down into five steps. The first line (1) represents the initial system before the iteration. We start with two system qubits (orange circles) and a bath (blue triangles). The first operation (line 1-2) is to add two new system qubits and two auxiliary qubit (red diamond). The new qubits are placed to the right of the old ones. The second operation (line 2-3) is to act with a unitary U between the indicated qubits and from (3) to (4) we apply another layer of unitaries and, thus, D = 2 for this example. Finally, in line (4)-(5) we discard the auxiliary systems. Evolution in the Heisenberg picture of an observable initially supported on the last two incoming qubits in the lower right corner (filled orange corner) at the fourth time step. The U indicates that we apply a unitary between the involved qubits and the dashed arrows indicate that old qubits are shifted to the left. Note that we suppressed the unitaries that do not contribute to the expectation value. The distance of the observables to the identity is indicated by how filled the element, i.e., empty shapes indicate that the observable is proportional to the identity on that system. Note that as in the Heisenberg picture we discard qubits, this causes the observables to mix. Moreover, we see that after two iterations the observables are essentially proportional to the identity and it suffices to implement that part of the circuit to estimate it. qubit in between every existing qubit and nearest neighbors interact, resulting in a tree structure.
Past causal cone Formally, the final state of the system can always be written as Here A t adds the new subsystems and auxiliary qubits, U t is a unitary channel that consists of D twoqubit gates for each edge in G t and D t traces out the ancillary systems. An important property of our framework is that it allows to bound the number of qubits that can influence the value of a local observable, referred to as the causal cone of an observable 18,19 . The growth of the casual cone depends on the geometry of the graph G t . To see this, consider an observable O T on the final state ρ. According to Eq. (1), the expectation value is Here Φ Ã t is the evolution in the Heisenberg picture.
We can use our framework to bound the size of the radius of the support of O T on the final state i.e. the number of qubits on the final state that O T involves. Going back to the t'th iteration, we denote the be the radius of the smallest ball in G t containing the support of O t . That is, O t differs from the identity on qubits that are at most 2R(O t ) edges away in the graph G t . To analyze the growth of the support and its past causal cone we consider the action of T acts by tensoring the identity operator on the auxiliary qubits, not increasing the support. In the next step, U Ã T increases the support. As O T has radius R(O T ), it will be mapped to an observable with radius at most R(O T ) + D by U Ã T according to the locality assumptions of U Ã T (i.e. the restriction of D two-qubit gates for each edge). The map A Ã T will then map this observable to O T−1 supported on qubits that correspond to vertices in G T−1 , as it traces out all the qubits added at iteration T. This can potentially decrease the support of the observable, as in DMERA.
Given the graphs G t and a constant D, it is straightforward to track the support of the observable and the past causal cone with the above procedure. This allows us to find the maximum number of unitaries (N U (t, r)) and qubits (N Q (t, r)) in the past causal cone of an observable with radius of support of r on the final state going back to iteration t. Note that N Q (t, r) keeps track of the total number of qubits necessary to implement the past causal cone and thus also includes those that were discarded at a previous step.
Estimating local observables So far we have devised a way of keeping track of the unitaries in the past causal cone of local observables. However, we are also interested in quantifying how much each iteration of the past causal cone contributes to the expectation value. In case the expectation value of the observable stabilizes after a couple of iterations, we can find smaller quantum circuits than the entire causal cone that will approximate the desired expectation value.
Inspired by refs. 12 Here ∥⋅∥ ∞ is the operator norm. The mixing rate, δ(t, r), quantifies how close observables on the final state, whose support is contained in a ball of radius r, are to the identity after going back to the t'th iteration of the evolution in the Heisenberg picture (Fig. 2). Intuitively speaking, δ(t, r) measures how many steps of the circuit contribute to the expectation value of local observables before it stabilizes, as we are interested in the regime in which c approaches tr ρO ð Þ for large enough t. This is also connected to the memory of the evolution 24 . The next lemma formalizes this intuition (see "Methods" for a proof): Lemma Let O T be an observable supported in a ball of radius r. Then , which holds for all ρ 0 . In other words, only the last T − t steps of the circuit are necessary to approximately compute the expectation value of O T up to an error of 2δ(t, r)∥O T ∥ ∞ . Note that the expectation value is independent of the initial state ρ 0 , which we furthermore may restrict to the qubits that are in the support of O t . We may further reduce the size of the circuit that needs to be implemented by restricting to the past causal cone. Combining these two observations leads to the statement of our main result: Theorem Let O T be an observable supported in a ball of radius r and ρ 0 be a state on the qubits that are in the support of O t . It is possible to compute tr ρO T ð Þ up to an additive error 2δ(t, r) by implementing a circuit consisting of N U (t, r) two-qubit gates on N Q (t, r) qubits.
The theorem implies a way of performing VQE given bounds on δ(t, r) using a potentially smaller NISQ device than preparing the whole state. This is because implementing the smaller, effective circuit, requires fewer qubits and gates.

Robustness to noise
Consider the objective of calculating the ground state energy of a two local Hamiltonian H, in the sense that it only acts on nearest neighbors in G T , with each local term H i satisfying ∥H i ∥ ∞ ≤ 1. It suffices to estimate all H i individually to obtain an estimate of the global energy of the state by adding up the energy terms. Now suppose that we can implement each 2-qubit gate with an error ϵ U in operator norm and can prepare each initial qubit up to an error ϵ P . This implies that the total error of implementing the causal cone and measuring each H i is bounded by ϵ U N U (t, 2) + ϵ P N Q (t, 2). Thus, by only implementing the circuit from iteration t to T, it is possible to estimate the energy of each term with an error of 2δðt; 2Þ þ ϵ U N U ðt; 2Þ þ ϵ P N Q ðt; 2Þ: This generalizes to observables with arbitrary radius r and can be improved to by exploiting the fact that : 2δðt; rÞ þ P T k¼tþ1 δðk À 1; rÞϵ U N U ðk; rÞ À N U ðk À 1; rÞ ð Þ þ P T k¼tþ1 δðk À 1; rÞϵ P N Q ðk; rÞ À N Q ðk À 1; rÞ ð Þ : To see this, recall that δ(k, r) measures how close the operator O k is to being proportional to the identity, as there exists an operator A k with the same support as O k such that O k can be decomposed into O k = c1 + δ(k, r)A k . At the k'th iteration, any evolution in the Heisenberg picture only acts non-trivially on A k and changes the expectation value of the observable w.r.t. to any state by at most δ(k, r)∥O T ∥. Thus, if we actually implement a noisy version of the original evolution which is ϵ U close to it, then we can only notice the effect of the noise in the part given by δ(k, r)A k . We conclude that each noisy unitary contributes with an error at most ϵ U δ(t, r), i.e., the effect of noise decreases in time if δ(t, r) decays. As there are N U (k, r) − N U (k − 1, r) new unitaries in the causal cone at the iteration, we obtain the bound. We note that the noise robustness we obtain is of the same order as the one obtained implementing the whole circuit, as in refs. 12,13 , and refer to Supplementary Information for a detailed discussion. These results are related with the fact that δ(t, r) and the geometry of the interactions govern the correlations present in the state produced. For E T , F T two observables of disjoint support of radius r and t be the largest t such that E t and F t have supports that intersect we can show that: tr ρE T F T ð ÞÀtr ρE T ð Þtr ρF T ð Þ j j 6δðt; rÞ: As the decay of δ(t, r) also governs the noise robustness, we see that there is a trade-off between the correlation length and the robustness to noise. For instance, one should expect δ(t, r) to be exponentially decaying for states with finite correlation length.

Estimating convergence
It is necessary to bound δ(t, r) in our approach in order to bound the error. Thus, it is important to find conditions that guarantee the decay of the mixing rate and to develop protocols to estimate the mixing rate on a NISQ device. In the translationally invariant case, one can apply the large toolbox available to estimate mixing time bounds [25][26][27][28][29] , as further explained in Supplementary Information.
However, it is important to acknowledge that obtaining rigorous mixing time bounds is notoriously difficult even for classical systems 30 . But this has not kept Markov Chain Monte Carlo algorithms from being one of the most successful methods to simulate physical systems 31 . There exist many heuristic methods for classical systems 32 and here we also discuss a heuristic method to determine when the circuit has stabilized. As made transparent by Eq. (2), whenever the circuit has converged, the output of the circuit is independent of the initial state ρ 0 . Thus, one possible way of checking that the circuit has indeed converged is testing several different initial states and making sure that the expectation value of the output with respect to the observable does not depend on the input. This can be done by picking a set of initial states that is overcomplete, i.e., spans the space of all states as detailed in the "Methods" section.
In short, the approach is to pick random initial product states on the support of the observable O t and compare the value of the expectation value to that of the initial state where all qubits in the support are in state 0 j i. If the expectation value with respect to different initial states all coincide, we build confidence that the computation has indeed converged. On the other hand, if the expectation values differ for two different initial states, then we have not converged and must go deeper and decrease t. We denote the maximal difference of the expectation value for the several different choices of initial state by Δ z and refer to the "Methods" section for a precise definition. An example of the approach is shown in Fig. 3 for a matrix product states with both fast and slow convergence.
As we can see from the figures, estimating Δ z gives reliable convergence diagnostics. Importantly, this is obtained with only a modest number of randomly selected input states. This suggests that estimating the convergence does not outweigh the overall advantage of bounding the circuit size for estimating local observables with our framework. We do note, however, that this is a heuristic and not rigorous approach. For rigorously establishing convergence, it would in general require a sample complexity increasing exponentially with the support. We refer to the supplementary information for more details.

DISCUSSION
To demonstrate the implications of our results, we summarize the noise robustness and required number of gates and qubits in Table 1 for some interaction schemes. We are able to significantly decrease the number of unitaries and qubits compared to the approach of refs. 12,13 . This is because we only require the circuit corresponding to the past causal cone until it stabilizes to be Fig. 3 Convergence check. Illustration of convergence check for MPS in a rapid and slow mixing scenario. We consider a bath of 9 qubits arranged on a line. In each iteration a new system qubit is added and a fixed circuit is run between the bath and the new system qubit. In the rapid mixing case, we perform a depth 3 circuit between the system qubit and the bath such that each gate is followed by a depolarizing channel with probability 5%. The complete state was evolved for 50 time steps and, thus, consists of total of 50 qubits, and we measured Pauli observables on the last two qubits. In the slow mixing we perform a depth 2 circuit instead. If we denote by X t the expectation value of the state evolved by t steps, we define the relative error to the true value to be |1 − X t / X 50 |, as X 50 corresponds to preparing the whole state. To generate each plot, we have generated 50 instances where the gates in the circuits were picked at random for each instance. The figures display the logarithm in base 10 of the average of the quantities. Even with moderate levels of noise, we can faithfully reproduce the expectation value up to a relative error of 10 −2 after 5 iterations, giving an order of magnitude saving in the number of qubits in the rapid mixing scenario, while in the slow mixing, we required around 13. Moreover, we have picked 20 random initial states to estimate Δ z , as more initial states did not seem to improve the estimate. As we can see, Δ z is a good proxy for the relative distance to the true expectation value.
Error and resources required for implementing the effective causal cone, as in (2), for different interaction schemes under the assumption that δ(t, r) = ce −λ(T−t) . RI-d refers to the case in which a d-dimensional bath interacts with a d-dimensional system at each iteration. All entries are only up to leading order in D, T and λ. For λ independent of system size, we see that it is possible to approximate all expectation values with quantum circuits of constant size.
implemented, in contrast to the the whole causal cone. Clearly, these results also imply that it is possible to approximate these expectation values classically if the resulting effective circuits are of a classically simulatable size. Our results provide an intuitive understanding of the stability of these computations. Each iteration contributes less to the value of expectation values, which implies that there is a small effective quantum circuit underlying the computation. Furthermore, the size of this circuit is related to the correlation length of the state and the effect of noise decreases proportionally to the correlations between regions.
Although rigorously testing at which iteration the circuit has converged might require exponential resources, we see that chosing a few random initial states and comparing the different expectation values provides useful guidance to check wether convergence has occurred. This shows that it is feasible to build confidence in the convergence and required depth with moderate resources using such heuristics. However, it would still be interesting to establish more rigorous protocols under suitable additional assumptions.
There are significant challenges in scaling up current qubit technologies [33][34][35] . The reduction in the number of qubits that we have shown above means that it may be possible to explore many-body quantum states with NISQ devices with substantially fewer qubits, potentially bringing such tasks into reach for current technology 3,4 . The possible reduction in the number of gates also reduces the necessary runtime of the circuits, which is important for hardware subject to qubit loss over time such as trapped atoms 36 . Note that the objective of this method is to expand the simulating capabilities of NISQ devices subject to strict hardware limitations. This is in contrast to other techniques, like measurement regrouping 8,[37][38][39][40][41] , that focus on optimizing resources given the ability to implement the whole circuit that prepares a given state.

Proofs and checking mixing
The main result of our work is based on the lemma in the main text. In order to prove this lemma, we have used a method based on viewing the generation as a quantum channel in the Heisenberg picture. The formal proof is given below.
Proof By the definition of δ(t, r), we see that where A is some observable supported on supp(O t ) satisfying ∥A∥ ∞ ≤ ∥O T ∥ ∞ . As Φ Ã ½0;tÀ1 is a quantum channel in the Heisenberg picture, Φ Ã ½0;tÀ1 ð1Þ ¼ 1 and k Φ Ã ½0;tÀ1 k 1!1 1 42 . Thus, where k Φ Ã ½0;tÀ1 A ð Þk 1 k O T k 1 . We conclude that 2δðt; rÞ k O T k 1 by combining (5) and (6). □ Our results show that having an estimate for the mixing rate δ(t, r) allows to bound the number of qubits and circuit size needed to estimate local observable on a many-body quantum state. While certain classes of states are known to have rapid mixing leading to fast convergence of the observables (like ground states of gapped Hamiltonians), we have developed a heuristic method for estimating the convergence of local observables. The method relies on the observation that we can expand any density matrix using a basis of single qubit states. Let ρ t = Φ [0, t] (ρ 0 ⊗ ρ b ) be the state we obtain by evolving from time 0 to t and O T be defined as usual. Moreover, let m = N Q (r, t) be the number of qubits in the support of O t . To check the convergence of the circuit, we prepare input states ψ z j i¼ N m i¼1 z i j i, where z i j i 2 f 0 j i; 1 j i; þ j i; À j ig. The states ψ x j i thus corresponds to various product state combinations of the states 0 j i; 1 j i; þ j i and À j i on the support of the observable O t . It is wellknown that these states form a basis of the space of Hermitian matrices. We can therefore expand ρ t as where a z 2 C. Furthermore, it is easy to see that ∑ z a z = 1 by taking the trace. Now define Δ z to be given by: From this, we immediately obtain that: Equation (7) suggests the simple protocol of picking random initial product states ψ z j i and comparing the expectation value with the outcome for initial state 0 j i m to estimate the convergence. If the expectation value with respect to different initial states ψ z j i all coincide with the one of 0 j i m , we build confidence that all Δ z are small, thus ensuring that the expectation value is similar as picking the initial state to be 0 j i m .

DATA AVAILABILITY
No datasets were generated or analyzed during the current study.