Preparation and verification of tensor network states

We consider a family of tensor network states defined on regular lattices that come with a natural definition of an adiabatic path to prepare them. This family comprises relevant classes of states, such as injective Matrix Product and Projected Entangled-Pair States, and some corresponding to classical spin models. We show how uniform lower bounds to the gap of the parent Hamiltonian along the adiabatic trajectory can be efficiently computed using semi-definite programming. This allows one to check whether the adiabatic preparation can be performed efficiently with a scalable effort. We also derive a set of observables whose expectation values can be easily determined and that form a complete set, in the sense that they uniquely characterize the state. We identify a subset of those observables which can be efficiently computed if one has access to the quantum state and local measurements, and analyze how they can be used in verification procedures.


I. INTRODUCTION
The simulation of many-body states is one of the most promising and long-awaited applications of quantum computing. In particular, quantum computers are expected to efficiently prepare certain entangled multipartite states, which can help us in the study of quantum many-body systems, or in variational quantum algorithms. The advent of the first generations of both analog and quantum computers has triggered a strong interest in the preparation of such states. For instance, GHZ states up to tens of qubits have been prepared with trapped ions [1][2][3][4][5], Rydberg atoms [6], superconducting qubits [7][8][9][10], photons [11][12][13][14] or nuclear spins [15,16].
Tensor network states (TNS) constitute an especially appealing family of multipartite states [17]. On the one hand, they are expected to efficiently approximate the ground states of local Hamiltonians. On the other, many paradigmatic states in the realm of quantum information or condensed matter physics are simple examples of TNS. The best-known class of such states is that of Matrix Product States (MPS) [18], which corresponds to a one-dimensional geometry. Higher-dimensional generalizations are known as Projected Entangled Pair States (PEPS) [19]. In both cases, they are characterized by the bond dimension, D, which is directly related to the size of the tensors building such states. Those states possess a special property, namely that they are the ground state of a local, frustration-free Hamiltonian. This implies that, in case of no degeneracy, one can easily check the successful preparation of the state by measuring a set of local observables. Thus, such states naturally play an essential role in the certification of quantum computers [20][21][22][23][24][25][26][27][28][29][30][31][32].
The preparation of TNS has been actively pursued in the last years and, in particular, methods that operate efficiently in terms of the number of qudits (or lattice size), N , have been devised. 1 Matrix Product States can be sequentially generated [33] in a time that scales linearly with N . In fact, MPS of up to ten qubits have been recently prepared in a superconducting setup [34][35][36][37]. Certain kinds of PEPS (the so-called sequentially generated) can also be prepared in the same time scale [38] and proposals for the generation of sequentially generated PEPS have been recently put forward [39,40]. While containing many paradigmatic examples of TNS, those states are fine-tuned in the sense that a small change on the tensors defining the state may lead to another PEPS outside that class. In fact, those tensors are strongly restricted by the fact that the states have to be sequentially generated.
In [41], a very efficient quantum algorithm to generate a wide range of PEPS was introduced. That class of states is stable under deformations of the tensors and thus they are not fine-tuned. The algorithm is based on the adiabatic method and the circuit depth scales as O(log(N )). The algorithm needs, however, the existence of a gap, ∆, along the adiabatic path, something which is difficult to ensure since checking that typically requires a computational time that scales exponentially with N . Additionally, it is devised for digital quantum computers, but not analog ones.
In this paper we consider a family of states on arbitrary lattices and prove two novel results on this family. First, we show how the computation of the gap of the parent Hamiltonian can be expressed as a semidefinite programming problem (SDP), and how this allows us to efficiently compute lower bounds δ ≤ ∆ on the gap. Second, we show that for such families of states, it is possible to predict the expectation values of many observables beyond those appearing as terms of the parent Hamiltonian.
In fact, there is an exponential number of such observables, which forms a complete set in the set of operators acting on the many-body Hilbert space.
The first result naturally leads to preparation protocols by using for instance the one in Ref. [41], since we can efficiently identify the subset of tensor network states for which a bound in the gap is known. Besides, we also extend the adiabatic algorithm to continuous time, which is more suitable for analog quantum computers.
The second result naturally leads to certification protocols based on interactive proofs (see [20][21][22][23][24][25][26][29][30][31][32] for previous works on interactive verification schemes) which are inspired by the difficulty of sampling from PEPS [42,43]. While they require a certification time that grows exponentially with N [31,32], we propose efficient versions that, however, rely on stronger standard complexity assumptions. We also explain that, in case they can be spoofed, this would immediately lead to new classical algorithms to estimate physically relevant expectation values of observables in the class of states we consider.
This paper is structured as follows. In Section II we present the class of states and their corresponding parent Hamiltonian. The states depend on two positive parameters, t, β ≥ 0, that can be viewed as time and inverse temperature, respectively. We show that, by construction, these states can be smoothly connected to a product state. This family of states is esentially the same as the one considered in Ref. [41], although our formulation allows to extend the adiabatic quantum algorithm of [41] to continuous time in a natural fashion. We will see later how this family includes injective MPS and PEPS, as well as some ground states of classical models whose PEPS description might not be injective (although their parent Hamiltonian has unique ground state). In Section III we first show how one can efficiently find lower bounds on the gap of the parent Hamiltonian by means of an SDP. In particular, for every value of t we can find a maximum value of β(t) such that for all β < β(t) the gap of the Hamiltonian can be lower bounded by a constant that does not depend on the system size. We then introduce sets of operators whose expectation values can be easily computed and that are complete, in the sense that they provide a tomography of the state. Finally, we propose verification protocols in Section IV and discuss some possible complexity arguments.

II. STATES AND PARENT HAMILTONIAN
In this section, we introduce the family of states that we will consider in the present work, which is comparable to the one considered in Ref. [41]. It is built in terms of sets of local commuting operators, together with two parameters t, β ≥ 0. For t = β = 0, they are product states, whereas otherwise they are entangled and can be efficiently expressed as PEPS. Based on that fact, we will explicitly construct a frustration-free parent Hamil-tonian, which will play an important role in the procedure to prepare the states. We then analyse how to prepare these family of states adiabatically and study the scaling of the computational time as a function of the system size and a lower bound on the gap. We build upon the work done in Ref. [41] and extend their adiabatic algorithm to continuous times. We review here the main idea of the algorithm from Ref. [41] and its runtime, and refer to the reader to the original paper for a more technical discussion. Note that while in Ref. [41] a constant gap was assumed, we will present later in Section III a method for lower-bounding such gap, which constitutes our main result regarding ground state preparation. This will immediately allow us to know which states we can efficiently prepare. Lastly, we show how the presented family of states includes many physically relevant examples of TN states.

A. Setting
We will consider a rather general setup, although we will give examples later for regular lattices. We consider N qudits, with Hilbert space H d = C d and computational basis {|0 , . . . , |d−1 }, located at the vertices V of a graph G(V, E) with edges E of bounded degree z (that is, the maximal number of edges starting from a given vertex). We define H = H ⊗N d , and denote the set of Hermitian operators acting on H by A. For any O ∈ A, we define its support λ(O) ⊂ V as the subset of vertices such that where tr i is the trace with respect to the qudit at vertex i.
The graph G defines a natural distance d(i, j) between two vertices i, j ∈ V as the minimal number of edges connecting them. We also define the radius of a subset of vertices λ ⊂ V, We denote by A r ∈ A the set of Hermitian operators acting on H whose support has a radius of at most r.

B. States
We call K r,M ⊂ A the set of operators that can be written as where κ n ∈ A r , κ n ∞ ≤ 1 and κ n 0. That is, it is the set that can be written as a sum of M commuting, positive semidefinite, subnormalized operators that have a radius of at most r. Apart from the trivial cases, where the operators κ n act on single vertices or when they are products of the same single-qudit operator, one can easily construct non-trivial sets K r,M . In Appendix A we briefly review some of them.
Given r α , M α ∈ N, K α ∈ K rα,Mα , for α = 1, 2, we define the family of states where β, t ≥ 0, Z is a normalization constant and |ϕ i are arbitrary single-qudit states. This family of states obviously contains all product states if we take β = t = 0. For β = 0, we have Z(0, t) = 1. Note that we only explicitly denote the dependence of |Ψ(β, t) on β and t, while omitting the dependence on K α , for α = 1, 2, and |ϕ i , in order to ease the notation.

C. Parent Hamiltonian
We now show that any state (3) is the unique ground states of a local, frustration-free Hamiltonian, which we construct explicitly. We denote by κ α,n the operators appearing in the decomposition (2) of K α , for α = 1, 2. For j ∈ V, let µ j := {n | j ∈ λ(κ 2,n )} be the index set of all terms κ 2,n which act non-trivially on site j, and ν j = {m | ∃ n ∈ µ j : λ(κ 1,m ) ∩ λ(κ 2,n ) = ∅} the index set of all terms κ 1,m which overlap with one of the previous terms. See Figure 1 for an example of how such set of vertices would be. Then, for j = 1, . . . , N , we define where Π j = 1 j − |ϕ j ϕ j | j acts on the qudit at vertex j, and which is invertible. With this definition, we have h j = h † j 0 and h j |Ψ(β, t) = 0. We define the parent Hamiltonian of |Ψ(β, t) as Note that we have now also suppressed the dependence of h j on β and t for convenience. Let us show that, indeed, (3) is the unique ground state of such an operator. Since h j |Ψ(β, t) = 0, we have that H(β, t)|Ψ(β, t) = 0, and since H(β, t) 0, this implies that |Ψ(β, t) is a ground state of H(β, t), with ground state energy 0. Conversely, if H(β, t)|Ψ = 0, then h j |Ψ = 0 and thus Π j e −itK2 e −βK1 |Ψ = 0 for all j, which in turn means that |Ψ is proportional to |Ψ(β, t) : We thus see that |Ψ(β, t) is the unique ground state of H(β, t).

D. Adiabatic preparation
The existence of a smooth path of Hamiltonians connecting H(β, t), and thus |Ψ(β, t) , to a simple product state at H(0, 0) implies that these states can be prepared adiabatically, by starting with the product state |Ψ(0, 0) = |ϕ 1 ⊗ · · · ⊗ |ϕ N and adiabatically changing the Hamiltonian from H(0, 0) to H(β, t). In fact, the first step of the procedure -changing the Hamiltonian from H(0, 0) to H(0, t) -corresponds to applying a unitary transformation U = e itK2 to |Ψ(0, 0) . This transformation can be implemented exactly in time t by evolving with K 2 (rather than H). Alternatively, using the fact that K 2 is a sum of local commuting terms, U can be decomposed into a finite-depth local unitary circuit (where the number of layers only depends on the structure of the interaction), which can be realized exactly on a digital quantum computer or simulator in constant time.
The task that needs to be implemented adiabatically is the second part of the preparation, that is, the interpolation from |Ψ(0, t) to |Ψ(β, t) . Generally, the time required for a faithful adiabatic evolution will depend on the magnitude of the spectral gap of H(β , t) along the interpolation β ∈ [0; β]. As it turns out, for the given type of interpolation, we can devise an efficient way of checking the presence of such a gap numerically, which we present in Section III A. Once we have established a lower bound on such a gap, we can use any standard bound for adiabatic evolutions [44,45], which gives an adiabatic runtime scaling as T = O(N 2 ∆ −3 −1 ), with ∆ a lower bound on the gap, and the error in the final state.
Moreover, for the situation at hand, we can improve the scaling of the adiabatic preparation by making use of the locality of the Hamiltonian, combined with the version of the adiabatic theorem proven in Ref. [41]. To this end, we construct an alternative interpolation from H(0, t) to H(β, t) where in each step, we only change one of the terms in the Hamiltonian; importantly, the method derived in Section III A to prove gaps still applies in that case. By changing the imaginary time in a suitably smooth way along this interpolation, we then obtain that the adiabatic time required per Hamiltonian term changed scales logarithmically with the desired accuracy. We can now concatenate the interpolation for all the individual Hamiltonian terms. However, we would still have a Hamiltonian acting on the whole lattice, which will translate in an overhead in the runtime if we want to devise a digital algorithm. To overcome this, we use the fact that changes performed on distant terms can be equally well carried out in parallel due to an effective light cone through Lieb-Robinson-bounds [41,46]. This means that at each step we only work with Hamiltonians supported on a region of size poly-logarithmically in the system size. Combining all these results leads to a preparation scheme for |Ψ(β, t) where the adiabatic time scales poly-logarithmically with the desired accuracy and the problem size N , and thus exponentially better than other known methods for preparing MPS and PEPS. For a more technical discussion on this preparation method we refer to the reader to [41]. The bounds for the continuous-time version of the algorithm can essentially be derived from the ones presented in the original work.

E. Connection to tensor networks
Let us now show that the states (3) can be efficiently described as PEPS [17] with the same connectivity as G(V, E) and a bond dimension D which is upper bounded by a function of r 1 , r 2 , z, and d, but which does not depend on N or M . To see this, let us consider an edge (i, j) ∈ E. Since the operators κ 1,n commute pairwise, we can express the product of operators that act on both i and j as a single one, i,j∈λ(κ1,n) e itκ1,n = e it i,j∈λ(κ 1,n ) κ1,n . We can bound the number of terms that appear in the sum in the exponent by the number of operators that act on qudit i. Since the individual operators κ 1,n act on a radius r 1 , note that e it i∈λ(κ 1,n ) κ1,n acts on, at most, all the qudits that are at a distance less or equal than 2r 1 from qubit i. The number of such neighbors can be bounded by Finally, note that an operator that acts on x qudits increase the bond dimension at most by O(d x ). 2 Iterating this for every edge (which overcounts interactions) gives an upper bound O(d z 2r 1 ) for the required bond dimension. Since a similar argument is valid for the operators κ 2,m , we conclude that the states of the form (3) can be described by a tensor network of bond dimension at most O(d z 2(r 1 +r 2 ) ). Importantly, the bond dimension remains bounded independent of the system size since r 1 , r 2 and z are size-independent.
In the following, we particularize the above results to standard PEPS. In particular, we will show that it contains all injective MPS and PEPS, as well as a PEPS corresponding to classical models, which (e.g. on the square lattice) are described by non-injective PEPS with a unique ground state. On the other hand, we generally cannot expect that G-injective PEPS -that is, those which can exhibit topological order -are of this form, regardless of boundary conditions and thus ground space degeneracy, as they are not connected adiabatically to product states; however, this can be easily remedied by choosing an RG fixed point in the corresponding phase as an initial state, rather than a product state, allowing for all findings in this paper to be carried over with only minor adaptions.
We will in the following consider regular lattices in one or higher dimensions. In the first case, we have MPS, and for the higher dimensional case, we have PEPS.

Injective MPS
Matrix Product States (MPS) are the simplest TN [17,18]. They can be written as We consider a special subclass, so called injective MPS. They fulfill d = d 2 0 , and are constructed with two qudits L n and R n each of dimension d 0 per vertex, and FIG. 2: MPS construction in terms of the operator e iκ2,nt , that creates an entangled pair between sites R n and L n+1 , and e βκ1,n , that maps the virtual sites L n , R n to the physical qudit at position n.s where |Φ is a state where, at each vertex, the qudit R n is in a maximally entangled state with the qudit L n+1 on the vertex to its right (and thus L n in a maximally entangled state with R n−1 ), and 0 < Q n ≤ 1 1 are invertible operators that transform C d0 ⊗ C d0 → C d . Generically, MPS become injective by blocking ≤ D 4 original qudits [47]. MPS obviously fall within the family of states defined in (3) for the special case of a 1D lattice as graph. Here, for each n, κ 2,n acts on R n and L n+1 in such a way that e iκ2,n t creates the maximally entangled states on those qudits. The operator κ 1,n are given by κ 1,n = − log(Q n )/β and act on a single vertex each. This is illustrated in Figure 2.

Injective PEPS
Injective PEPS are the generalization of (9) to higher spatial dimensions. The state has the same form, but now there are z n qudits at site n, where z n is the coordination number of vertex n (i.e., the number of edges connected to the vertex). The state |Φ contains entangled pairs along all possible vertices of the lattice. As for MPS, one can readily see that these states lie in the family defined in (3).

PEPS corresponding to classical models
Let us consider a classical model in a lattice. These ground states can be described by PEPS, which might be non-injective, but whose parent Hamiltonian has a unique ground state [48]. Consider a Gibbs state of some classical Hamiltonian H cl (s 1 , ..., s n ) = where Z is the classical partition function. We can rewrite this state as follows [49]: whereĤ cl is an operator with eigenstates |s 1 , ..., s n and eigenvalues H cl (s 1 , ..., s n ). This state is of the type of Eq. (3), where K 1 = H cl /2 and K 2 = 0. A description of these states in terms of PEPS can be easily obtained.
For the special case of two-body nearest-neighbor interactions, we get a bond dimension equal to the dimension of the physical degree of freedom, see [49].

III. GAPS AND EXPECTATION VALUES
In this section we establish our two key results, that find application in both efficient preparation protocols for the states (3) as well as for their verification. First, we develop an effective method to compute lower bounds to the gap of the parent Hamiltonians (6) for some range of parameters t, β. This will ensure that the corresponding ground states can be efficiently prepared with the adiabatic algorithm presented in the next section. Second, we determine the expectation values of a complete set of operators in those states, so that one can use them to check that the state has been successfully prepared.

A. Gaps
We will now explain how to efficiently obtain a lower bound on the gap of the parent Hamiltonians constructed in the previous subsection, for specific points in parameter space as well as uniform bounds for a whole parameter regime. To start with, note that H(β = 0, t) is a sum of commuting projectors h j for any t, and thus has a gap ∆ = 1. It thus remains to supply methods to bound the gap in the case where β > 0.
Let us first discuss how to obtain such a bound for a specific point H(β, t) = h j in parameter space. To this end, consider the semidefinite program (SDP): For any feasible point of the SDP (and in particular the optimum in (12a)), summing (12b) over all i = j yields (counting each pair twice, and up to a factor 2) or equivalently which says that H has no eigenvalues in the interval (0; x). Since the ground state energy is 0 by construction, this implies the existence of a spectral gap ∆ ≥ δ above the unique ground state.
Since (12) is an SDP with the dimension of the constraints independent of N , it can be solved efficiently. Note that the SDP can be simplified considerably by setting and observing that (12b) is trivially satisfied in those cases, leaving a number of equations linear in N . Another possible simplification amounts to first solve for each in- i , so we can add positive contributions to both a ij and c ij while still satisfying (12b)); a bound on the gap can either be computed directly from y i or through the SDP (13) while keeping a ij fixed. Finally, note that for β = 0 (where the h i are commuting projectors), condition (12b) holds for any choice of a ij = c ij , and thus indeed gives δ = 1.
Better bounds on the gap can be obtained by relaxing (12b) to only hold when summed over specific groups of terms (which still implies (13)); a natural such case would be the variant of the SDP obtained by first grouping adjacent termsh i = h i (with the sum over terms in some neighborhood) and then setting up the SDP for theh i . Furthermore, one can replace h i (also after blocking) by projectors with the same range as h i , since those relatively bound h i and thus have a system-size independent gap if and only if the original Hamiltonian does.
Having described a way how to efficiently obtain a lower bound to the gap for a given point H(β, t), how can we use this to build methods for certifying a gap over a whole range of parameters? The idea is based on the continuity of the SDP conditions (which should come as no surprise, given the finite dimension and smooth dependency of β of all objects involved). Given a certified gap δ for some H(β, t) using the SDP (12) (with corresponding optimal parameters a * ij and c * ij ), we show in Appendix B that the SDP for H(β + τ, t) has a feasible point with a ij = a * ij , and where for 0 ≤ a ij ≤ 1, and a variant thereof if a ij < 0 or a ij > 1, see Eqs. (B13). Importantly, c ij changes uniformly continuous as τ is increased starting from τ = 0: Thus, this allows one to obtain a lower bound δ(τ ) on the gap of H(β + τ, t) by virtue of (12d); importantly, the lower bound is uniform for a given interaction geometry and independent of the system size or the specific chosen model, and changes uniformly continuous with τ . We now take the τ 0 for which the lower bound closes, i.e. δ(τ 0 ) = 0, and re-run the SDP for H(β + τ 0 , t): If that SDP returns a gap as well, this proves that H(β + τ, t) is gapped for all 0 ≤ τ ≤ τ 0 . By starting from β = 0 (for which (12) trivially holds), we can thus establish the existence of a gap (and obtain explicit lower bounds to it) for a finite regime 0 ≤ β ≤ β 0 by evaluating the SDP (12) only at a finite number of points (where the bound can be improved by increasing the density of points). Note that for translational invariant systems, we can show that the SDP from (12) can be "symmetrized" (see Appendix C). This means that for a given feasible point, with values of the variables (a * ij , c * ij ), we can find another feasible point by averaging a * ij , c * ij , such that all coefficients a ij and c ij are equal for all pairs of terms in (12a). Moreover, this solution provides the same value of δ. This directly shows that in this case δ does not change with the system size, thus recovering the result first showed in Ref. [18] about the existence of a gap independent of the system size. Therefore, for TI systems the size of the neighbourhood τ is independent of the system size.
Let us now briefly discuss the suitability of the method for the TN of the previous section. For translational invariant (or suitably uniform) injective MPS, it was proven by Fannes, Nachtergaele, and Werner [18] that the parent Hamiltonian is always gapped. Indeed, they showed that by blocking a number of sites proportional to the correlation length (for a fixed bond dimension), which in turn can be bounded as a function of the gap, and replacing the blocked Hamiltonianh i by projectors, Eq. (12b) is fulfilled (with all a ij and c ij equal, as mentioned in the previous paragraph) -and thus the SDP will yield a gap δ > 0 -already with the restriction (15) imposed. In higher dimensions, an analogous result has been obtained for PEPS which are unique ground states of local Hamiltonian (in particular, injective PEPS) [50]: Whenever the Hamiltonian is gapped in the thermodynamic limit, the SDP with condition (13) will be satisfied by the projector-valued Hamiltonian obtained after blocking a number of sites which only depends on the gap and the geometry of the system, but not on the system size or details of the model.
Finally, regarding the PEPS corresponding to classical models in dimensions higher than one, for sufficiently high temperatures (small β) the SDP method will give a gap (as the Hamiltonian at β = ∞ is trivial). Note that the corresponding classical model may have a phase transition at sufficiently low temperatures, which implies that the correlation length will diverges and thus the gap will vanish in the thermodynamic limit (N → ∞). Thus, the SDP automatically also allows to determine an upper bound to that critical temperature, which will yet again get more accurate as we consider larger regions, both by blocking and by relaxing (or omitting) the restriction (13). Thus, the continuity bound of Appendix B at the same time provides a means of determining upper bounds on the critical temperature of classical statistical models.

B. Expectation values
Computing expectation values of ground states of local Hamiltonians is hard in general [51]. However, for ground states of frustration-free Hamiltonians, there are certain observables for which it is straightforward to compute such values. In this section, we find a complete set of operators for which this can be done and which will be on the basis of the verification protocols presented below.
As argued in Section II C, for |Ψ given in (3) we have h n |Ψ = 0, and thus trivially Ψ|h n |Ψ = 0. We will now define a set of observables for which one can also compute the expectation value in |Ψ . We will restrict to qubits (d = 2) and will take |ϕ j = |0 , although it is straightforward to extend the method to qudits and other states. We will denote the Pauli operators acting on the j-th qubit, j = 1, . . . , N , by σ j α with α = x, y, z; The key idea is to notice that for any λ ⊂ V, where with µ(λ) = j∈λ µ j , ν(λ) = j∈λ ν j , and where the sets µ j , ν j have been defined in Section II C. The first set of operators is defined as -that is, |Ψ is a right (left) eigenvector of Z λ (Z † λ ), using (17). We then have and thus Z ± λ have fixed expectation values. The second set is more ample. Given any P supported in λ ⊂ V with the property that there exists j ∈ λ such that 0|P |0 j = 0, we define Again, using (17) we can compute the expectation value of those operators The set of operators defined above, taken jointly for all λ ∈ V, is complete in A in the sense that their expectation values completely determine the state. To show that, we just have to devise a subset of 4 N linearly independent operators from that set.
To this end, let us start from the set of all operators P which are a product of Pauli and identity operators, and associate to each of them an operator Q ≡ Q(P ) using one of the constructions above. Each of these operators can be written as where α j = x, y, z, and λ ⊆ V is the set of sites on which P acts non-trivially. In case α j = z for all j ∈ λ, we define Q(P ) = Z + λ ; otherwise, Q(P ) = Q λ , Eq. (21). Finally, for P = 1 1, let Q(P ) = 1 1. In this way, starting from all products of Pauli operators and the identity, we have obtained a set of operators Q n , n = 1, . . . , 4 N . This set is linearly independent iff the matrix B n,m = tr(Q † n Q m )/2 N is not singular. Trivially, for β = 0, B n,m = δ n,m and thus not singular. Since the operators O λ used to define the map P → Q are analytic in β > 0, the determinant of B n,m will be analytic as well, and thus it can only vanish at countably many points, all of which are isolated. Thus, generically it will be non-zero and, in the possible measure-zero cases where it is can be circumvented by taking a closeby value of β.
The fact that the set of operators Z λ and Q λ is (over)complete means that any observable can be expanded as a linear combination thereof. If we are able to obtain the corresponding coefficients, then we will be able to compute the expectation value of all physical quantities. In practice, this will be difficult since, due to the fact that the operators O λ non-trivially overlap with each other, we will typically need an exponential number of terms in the expansion. Nonetheless, there may be a way of truncating that expansion, which would give lead to new algorithms in terms of tensor networks. Furthermore, since we know the expectation values of a basis of operators, we possess full tomographic information on state. However, as before, it is not useful to compute other expectation values. In Appendix D we show that the norms of the observables Q λ decay at most exponentially with their support size, and thus they are guaranteed to have bounded (polynomially-decaying) norm when the size of the support is fixed (at most O(log N )).
The operators Q λ are supported on a set of vertices that is larger than λ (roughly speaking, on all vertices that lie at a distance up to 2r from that set). It is possible to define other observables which have a smaller support and for which we can still compute the expectation values. This is relevant for more practical applications, like the verification protocols introduced below, where we want to make statements about measurements performed within the support of Q λ , and we want them to include as few qubits as possible. For that, given j ∈ V and an operator P supported on λ with j / ∈ λ, we define which only act on the joint support of O j and P . Again, using (17) we find the expectation value of those operators to be Ψ|Q j,α |Ψ = 0.
In particular, if we choose P to be an arbitrary Pauli product, then Q j,1 and Q j,2 can be used to replace the operators Q λ from Eq. (21) in our complete set of operators. At the same time, the operators Q j,1 and Q j,2 are still products of Pauli operators almost everywhere, except on the support of O j which has a fixed size. This means that each O j,1 , O j,2 can be estimated efficiently (i.e. with a number of (Pauli) measurements which only depends on the accuracy but not on the system size), and the only remaining operators in the complete set which cannot be estimated efficiently individually are the Z ± λ . In summary, we have defined a set of observables whose expectation values are either zero or one. We can choose a subset thereof where |λ| ≤ c, with c a constant independent of N . In such a case, since we know O λ we can efficiently write those observables as linear combinations of Pauli operators in the support of O λ (like (19)), or even smaller (like (24)). Note that the norms of these observables will also be efficient to compute in these cases, as we show in Appendix D.

IV. VERIFICATION SCHEMES
In this section we discuss different ways of exploiting the state preparation procedure for the state (3) as a verification scheme. First, we will analyse a quantum state verification method and show how to certify that the state has been created successfully by performing local measurements and using the fact that there exists a parent Hamiltonian that is both gapped and frustration free [52]. Then, we will consider the scenario of classical verification of quantum computation, where the goal is to make sure that someone else is in possession of a quantum computer solely through classical communication [21,22,[29][30][31][32]. We will restrict here to the case of qubits, although the arguments can be easily extended to qudits.

A. Quantum state verification
Unique ground states of frustration free Hamiltonians, like the ones we are dealing with here, can be trivially certified with local measurements. This is achieved by just performing local quantum tomographies to make sure that the expectation values h j = 0 for all h n defined in (4). Indeed, if this is the case, then H = 0 [cf. (6)], and since H ≥ 0 and has a unique ground state, this implies that the measurement must have been performed on that state.
In practice, since we can only perform a finite number of measurements, the estimates for h j will not be exactly zero; additionally, measurement errors will give rise to errors in those quantities as well. However, one can still estimate the success probability of the preparation in different ways. The most straightforward one is to relate the obtained expectation value of H = tr(Hρ) (with the corresponding estimated error) to the overlap of the state we have prepared, ρ, and the target ground state |Ψ . It is straightforward to show that the fidelity obeys where ∆ is the spectral gap of H, and δ the lower bound obtained in the previous section. Thus, if we can obtain this expectation value (with an error bar) and we compute δ, we directly obtain a lower bound to the overlap. Neglecting measurement errors, for a finite number of measurements, the estimate of Ψ|ρ|Ψ will have some error bar. To use the bound (26), we need that the error in tr(Hρ) = H = h j is sufficiently below δ. Assuming that we perform independent measurements and thus that we have independent errors with same standard deviation σ for the estimator of all the h j , we have that the error j for each term must satisfy δ/ √ N ∼ j ∼ σ/ L j , where L j is the number of measurements performed to estimate h j . If, in addition, we assume L j = L j , we obtain a conservative estimate for the total number of measurements of L tot N L j ∼ σ 2 N 2 /δ 2 . Note that instead of performing independent measurements, one could measure qubits belonging to non-overlapping regions in parallel [53], which might lead to a reduction of the total number of rounds.
The verification can also be analyzed as an adversarial game, where the prover prepares a state, ρ, gives it to verifier, and claims that it is indeed Ψ. The verifier can then perform measurements to gain confidence that this is indeed the case. Such a protocol as been analyzed in Ref. [28] by assuming that the verifier can measure by projecting in the basis of eigenstates of the local operators of the Hamiltonian h j . In Appendix E, we perform a similar analysis, but assuming that the verifier can only perform Pauli measurements, which might be a more realistic assumption for current experimental setups.

B. Classical verification of quantum computation
Let us now turn towards a different kind of verification, in which the verifier is fully classical and communicates with the quantum prover through a classical channel [21,22,[29][30][31][32]. Such protocols can also be formulated as an adversary game: Here, the prover claims that he can efficiently carry out quantum computations on a quantum computer. The verifier has to make sure that this is the case by communicating classically with the prover, that is, asking him to perform certain tasks on his quantum computer and report the results. Of particular interest is the case where both prover and verifier have limited additional resources, namely they can perform classical computations with a computational time that scales at most polynomially with the number of qubits. In this setting, the verifier can pose challenges to the prover which he can only accomplish if he has a quantum computer, but not with his limited classical resources, and the challenge is to find a way which allows the verifier with her limited classical resources to verify this.
Recently, several protocols achieving this task have been proposed whose security is based on standard complexity assumptions [29][30][31]. While the first protocol [29,30] is most adequate for fault-tolerant quantum computers, the latter [31] is very attractive since it can already be used to verify existing NISQ devices [54]. However, it requires the verifier to carry out classical computations whose runtime scales exponentially with the number of qubits, though with a relatively small exponent which makes it comparatively practical. Apart from their use to certify quantum computers that can be only used remotely by classical means, one of the most appealing applications of such protocols is in the context of certified random number generators [30].
Ground states of frustration-free quantum Hamiltonians that can be efficiently prepared, like the ones presented in this work, may offer an alternative way for this kind of verification; specifically, one can exploit the task of reproducing correctly the expectation values of the observables introduced in Section III B as a challenge for the prover. In the following, we will describe such verification protocols, and analyze their security and the underlying complexity theoretic assumptions in different regimes. As we will see, the straightforward application of this idea requires exponential time. More efficient versions of the protocol are possible, but the underlying complexity assumptions are less tangible and thus, their security remains unclear.

Protocol
The verification protocol consists of three steps: (i ) The verifier sends the prover instructions for preparing the state |Ψ . (ii ) This step consists of R rounds: in each round, the verifier sends the prover a set of ob-servables; the prover then prepares the state |Ψ , measures the observables, and reports the outcome. (iii ) The verifier performs certain checks on the accumulated measurement outcomes to verify that the prover has indeed prepared the state |Ψ and is thus in possession of a quantum computer.
For step (i ), the verifier just has to give the circuit that prepares the state to the prover, which she can e.g. obtain by trotterizing the adiabatic evolution. Alternatively, she can directly give the time dependent Hamiltonian H(t, β), together with the initial states |ϕ i , to the prover, e.g. in case he is in possession of an analog quantum computer. An honest prover will then be able to efficiently prepare the state |Ψ in a time that scales as log(N ), as discussed in Section II D.
For each round of step (ii ), the verifier sends the prover a list of bases α = (α 1 , . . . , α N ), α j = x, y, z, in which the individual qubits should be measured; the α will be generally drawn at random from some distribution which is dictated by the verification step (iii ). The prover then prepares the state and measures qubit j in the Pauli basis α j , and obtains results s = (s 1 , . . . , s N ), s j = ±1. For an honest prover, these results are drawn from a distribution After receiving the measurement basis, the prover prepares |Ψ and performs the measurements. Importantly, since |Ψ can be prepared in time O(log(N )), each round can in principle be carried out in time log(N ) as long as the prover has the ability to measure and communicate the results in parallel.
Step (ii ) allows for several natural generalizations. In particular, we can allow for measurements beyond the Pauli basis, we can enable the verifier to make adaptive queries, where the choice of s in any round can depend on the results of the previous round, and we can split each round into several sub-rounds of communication, where the state is prepared once and then a sequence of adaptive measurements is performed on it, measuring only a subset of the qubits at each time.
In step (iii ), the verifier uses the samples obtained from the prover to compute certain quantities which serve to verify that the prover is indeed sampling from the correct distribution P 0 . To this end, the verifier can e.g. use some of the quantities Q constructed in Section III B for which Q = 0, 1 (or variants thereof), or the Hamiltonian terms h j for which Ψ|h j |Ψ = 0 (which can be used to replace the Z − {j} ). Let us now consider one such observable Q supported on λ ⊂ V. It can be decomposed as where γ = (γ 1 , . . . , γ |λ| ), γ j = x, y, z and o(γ) are the expansion coefficients. The verifier can then estimate Q using an estimator wheres(γ) is the average outcome where the prover measured spin j in the basis σ γj for all j ∈ λ, and an arbitrary basis for all j / ∈ λ; that is, if we denote the set of rounds where α j = γ j for all j ∈ λ by R(γ), and the result of the r'th round by s r = (s r 1 , . . . , s r N ), we havē If the samples have been taken according to P 0 ,Ō → Ψ|O|Ψ in the limit |R(γ)| → ∞, which can be used as a means of verification (see Appendix E for a quantitative analysis). Note that alternatively, we can determine the average also using only rounds R(γ) where sites j / ∈ λ have only been measured in some specific bases.
As a simple example of this verification procedure, consider a resource state for measurement-based quantum computation (MBQC), such as the cluster state [55], or generalizations which allow for measurement-based computation using only Pauli measurements [56]. For those states, the preparation step is particularly easy, as they can be prepared by a single layer of commuting unitaries (i.e. β = 0) from a product state. They are frustruation free ground states of local Hamiltonians H = h i ≥ 0, h i |Ψ = 0, a property which allows for easy verification from the measurement statistics. And finally, we know that adaptive measurements in a suitable basis (Paulis and π/4 rotated states in the xy plane for the cluster state, or only Pauli measurements for the aforementioned generalization) are universal for quantum computation, making it impossible for a prover not in possession of a quantum computer to produce the correct distribution in such an interactive protocol, and at the same time imposing limitations on potential cheating strategies, as we will discuss in the following.

Complexity Analysis
Let us now analyze under which conditions the protocol allows the verifier to conclude that the prover must indeed be in possession of a quantum computer, what potential classical cheating strategies might be, and what limitations to such strategies exist.
a. Sampling from the quantum distribution is hard. The first cheating strategy would be to find a way to classically sample from the correct distribution P 0 , Eq. (27) -in that case, the verifier would have no way of detecting this. However, there are several strong complexitytheoretic arguments against that. First, note that -as already mentioned above -the adaptive version of the protocol contains measurement-based quantum computing: The cluster state is clearly in the given class (with β = 0), and an adaptive protocol with poly(N ) queries per round would implement a general quantum computation. Thus, sampling from the correct P 0 is impossible unless BQP = BPP. However, it is also known that sampling from the output of a circuit of commuting gates with random product states as input and measuring in a fixed basis (a case which is contained in our protocol for β = 0) is computationally hard, as shown in [57]. In order to prove the hardness of this setup, in [57] they relate the complexity of sampling from such a circuit with the complexity of sampling from IQP circuits, a task that is known to be hard unless the polynomial hierarchy would collapse to its third level [58,59]. We thus conclude that there is complexity-theoretic evidence that it is impossible for the prover to classically sample from the correct distribution P 0 (up to constant error) in polynomial time.
b. Reproducing all Q is hard. The second cheating strategy would be to sample from some different distribution P which is chosen such that all estimators for Q computed from P are correct. Any estimatorQ, Eq. (29), supported on a subset λ of sites can be computed in many different ways, which differ by the measurement settings over which we average for the sites not contained in λ, see Eq. (30) and the discussion below it. We will thus additionally impose that the estima-torQ converges to the same value for all those ways to compute it; this can be easily ensured by the verifier by computingQ in all different ways, or just in a randomly chosen one. (Alternatively, this can be ensured by computing the marginal probability distributions in different ways and checking their consistency, that is, by checking that the measurement setting α i of qubit i does not affect the distribution of the measurement results s j for any of the other qubits j = i; this can be seen as a kind of non-signalling condition on the distribution.) These conditions however imply that P = P 0 , the correct distribution derived from the quantum state. The reason is that for any given α, we can reconstruct P ( · |α) from the expectation values of all Pauli measurements given by arbitrary substrings of α (that is, where the Paulis in some positions have been replaced by identities, as in Eq. (23)) through an N -fold Hadamard transform (using the consistency condition above); the latter, in turn, can be reconstructed from all Q as they are related by an invertible transformation, as discussed in Section III B.
We thus find that this is not a viable cheating strategy, since no other distribution which yields all the correct Q exists. Note that together with the hardness results mentioned under point IV B 2 a, this implies that even sampling from a P which approximates Q for all Q sufficiently well is hard. The required accuracy in Q has to be chosen such that the expectation values of Pauli strings have exponential accuracy in N ; since the number of samples required to this end is in fact determined by the latter (the transformation (29) is exact), this generally requires an exponential number of samples.
Let us now discuss two different cases for the verifica-tion protocol and their security. First, we consider the case where both prover and verifier have bounded (polynomial) storage space, but the verification procedure can take an exponential number of rounds, and the verifier can use exponential time. In this case, the arguments from before regarding the hardness of sampling from P 0 or reproducing the Q imply the security of the protocol. Second, we consider the case in which the verifier only tests local observables, i.e., those supported on a region of at most O(log N ). Note that the verification is efficient in this case, i.e., it can be carried out in polynomial time for a prover and verifier having both polynomial storage space as before. In this case, we also argue how successful cheating strategies will be likely to fail. Most interesting, if such a cheating strategy succeeded, it would imply the existence of classical algorithms for computing local observables for generic tensor network states. c. A secure space-bounded and exponential-time verification protocol. In each round, the prover is only granted poly(N ) time (or logarithmic time with suitable parallel processing power). By performing an exponential number of queries, the verifier can get an exponentially precise estimateQ for any of the observables Q (either a randomly selected one, or all of them); importantly, the expansion coefficient o(γ) of Q in the Pauli basis, Eq. (28), is a trace of a product of local operators and can thus be computed in PSPACE, and the summands in (29) can be sampled and added sequentially one by one; the verifier thus only requires polynomial storage space. The prover, one the other hand, cannot classically sample from the correct distribution in the available polynomial (or logarithmic) time, due to the hardness results discussed in point IV B 2 a; at the same time, his limited memory prevents him from pre-computing all possible outcomes after having learned the adiabatic circuit for preparing |Ψ in step (i ) (even though the exponential time used for the overall protocol would allow for it).
Let us note that what makes our setup special is not the fact that the knowledge of all Q allows to reconstruct the probability distribution -such a complete tomography is possible in any scenario. Rather, what makes our construction special is that the verifier knows the required expectation values for all operators Q right away, whereas in a general tomography scheme, a computationally costly reconstruction procedure is required to obtain P from measured expectation values. However, in order to compute a general expectation value Q , the verifier still needs to collect an exponential amount of data from the prover: While it is sufficient to sample a small number of randomly chosen Q's, most Q's still have an exponential number of terms in their Pauli expansion (28). The situation here can thus be somehow regarded as the reverse from that introduced in Ref. [31]: While in that case, the verifier requires a polynomial number of measurements from the prover, it takes her an exponential time to check if they correspond to the correct probability (by computing the cross-entropy). In our case, the correct expectation values can be trivially computed, but one requires an exponential number of samples to obtain them. Note that computing the coefficients o(γ) of Q, which are needed to estimate its expectation value from the samples, might require exponential resources as well.
d. Proposal for an efficient verification protocol and its security implications. As we have seen, checking some Q , chosen at random from the set of all Q, allows to verify that the prover is in possession of a quantum computer. This, however, requires exponential resources, since a typical Q involves a significant fraction of the Pauli terms acting on its support, that is, an exponential number, and moreover we require an exponential accuracy. A practical verification protocol must thus resort to testing only local Q's, that is, those which are supported on at most order log(N ) sites, to an accuracy 1/poly(N ), as those can be sampled in poly(N ) time.
Thus, the question which arises is whether there is a way for the prover to efficiently sample from a distribution P which reproduces those local Q correctly, up to polynomial accuracy.
First, let us notice here that, if the prover indeed measures in the Pauli basis, reaching a polynomial accuracy on local Q expectation values uniquely identifies the ground state |Ψ among all multi-qubit states. The reason is that we can include the terms of the parent Hamiltonian H = h j in the set of Q's, and since |Ψ is the unique ground state of H and H is gapped, this implies |Φ − |Ψ < 1/poly(N ) if the h j are 1/poly(N )-close, see Eq. (26). Hence, to entirely establish the security of the protocol, it would be enough to prove that any distribution P reproducing Q with the desired accuracy must be obtainable by performing Pauli measurements on a multi-qubit state. Notice that a possible way to achieve that would be to prove a robust self-testing statement based on the Q expectation values (see [60] for a review of self-testing techniques).
Second, even if such a distribution P exists, there are constraints on the design of algorithms to sample from it efficiently. For instance, one might assume that it should be possible to characterize the space of solutions of the equations |Q− Q | ≤ 1/poly(N ), which locally constrain P , use them to find ways to sample locally correctly, and patch those ways of sampling together to obtain a sampling algorithm which works globally. However, such a strategy, if based only on the final conditions on Q , is bound to fail, unless further properties of P 0 are taken into account. Specifically, we can choose the Q to enforce 3-SAT clauses or some other classical NP-hard problem (such as spin glasses on a 2D lattice), in which case such an algorithm is bound to fail as it would have to solve the NP problem. (It is a fingerprint of hard instances of NP problems that the local constraints cannot be patched together easily.) Note that also knowing the precise local reduced density matrices does not necessarily make the problem easier, as the general quantum marginal problem is QMA-hard [61,62].
Thus, a successful cheating strategy most likely will have to use the full knowledge of the adiabatic path used to prepare the state, or equivalently knowledge of K 1 and K 2 . Note that Osborne showed that the computation of local expectation values on states that undergo an adiabatic evolution under a local gapped Hamiltonian can be achieved classically [63]. However, this method does not help the prover either, since it does not provide samples from the full probability distribution, which is what the prover is asked to return. Another approach could be to classically simulate the adiabatic evolution -for instance, one could try to adapt the Monte Carlo algorithm by Bravyi and Terhal [64] for the simulation of adiabatic quantum computation along a path which is both frustration free (which we have) and sign-problem free (which we don't necessarily have), followed by a measurement in the σ z basis (which we don't have). Indeed, it is plausible that such an algorithm will allow to correctly sample in cases where the final Hamiltonian is classical and only needs to be sampled in the σ z basis. On the other hand, it will likely break down if one of the conditions is not met, which will manifest itself in a sign problem in the Monte Carlo method. In fact, we cannot expect a cheating strategy based on simulating the adiabatic evolution to work, since such an algorithm (if working as desired) will pre-cisely sample from the correct distribution P 0 , which we have previously established to be computationally hard. Thus, in order to attack the protocol with such a strategy, one would have to devise a method to simulate the adiabatic evolution in a way where local expectation values are reproduced correctly, yet global properties would not (with the goal of circumventing hardness results); it is unclear how a route to accomplish this would look like.
take some Hermitian operators h A and h B acting on adjacent plaquettes A and B and assume that h A,B are of a Toric-code type [68] (see Figure 3a). By this we mean that = 0 for all α, α , β, β . Note that we can choose O(N ) free parameters, corresponding to the different evolution times {t α } α and {t β } β , as many as plaquettes in the lattice.
Bravy-type unitaries. Finally, we present a more general approach, which is a generalization of the results of Bravyi et al. [69] about the characterization of commuting, local Hamiltonians in 1D (see also [70] for a discussion on this topic). We will present the method for the particular case of a 2D rectangular locality structure, but the generalization to other geometries is straightforward. Let us start with a state defined over a lattice. At each site, we have four virtual particles of dimension D, which are mapped into a single physical qudit, with a physical degree of freedom d. For a visualization see Figure 3b, where we associate black dots with virtual particles; the big blue circles represent the physical sites, and correspond to the vertices of the lattice. We denote by P : (C d ) ⊗4 → C D the map from the virtual to the physical degrees of freedom. We take P to be unitary, which means that d 4 = D.
Let us denote by H q ∼ = C d the Hilbert space associated with the physical particle at vertex q. Bravy and Vyalyi [69] showed that H q can be decomposed as: where q ul , q ur , q dl , q dr correspond to the virtual degrees of freedom, see Fig. 3b. We denote by H q ul the Hilbert space associated with the virtual particle q ul and by H i q ul a subspace of H q ul . Let us now define some unitary operators U α , acting on four adjacent virtual particles, each belonging to a different vertex as in Fig. 3b. Here α just denotes an index to enumerate the different operators, without any spatial meaning. We do not impose any restriction on these operators apart from the virtual particles in which they act. Let us now show that when we consider the operators U α acting on the joint Hilbert space of two neighbouring physical sites, i.e. on H q ⊗ H q , they commute. To see this, consider the examples from Fig. 3b. We can write the operators U α and U α as: for some operators o qurq ul and o q dr q dl that act non-trivially only on the subspaces H qur ⊗ H q ul and H q dr ⊗ H q dl respectively. With this decomposition, it is easy to see that [U α , U α ] = 0. Now note that the application of the unitary operators P does not affect the commutation; it simply changes the basis in which the qudits are expressed. Hence we can define the final set of unitaries to be: where λ(U α ) denotes as before the support of the operator U α . Note that this defines a set of unitary operators over the full lattice, such that all of the operators commute pairwise. Note that this technique does not cover the generation of states where d = 2. Therefore for qubits other approaches like the ones presented above must be employed.
By construction, where (unless we have tighter bounds) α = e −2|νi\νj |∆τ , γ = e −2|νj \νi|∆τ , β = e −2|νi∩νj |∆τ (B5) (in particular, α, β, γ > 0). Henceforth, for convenience of notation we will write P := h i and Q := h j , as well as a := a ij , b := a ji , c := c ij , d = c ji . Let us now derive two inequalities. First, Second, using that 1 1 ≤ C −2 , −1 1 ≤ −γC −2 , −A 2 ≤ α1 1, and B 2 ≤ 1 1, we have that -for now assuming 0 ≤ a ≤ 1 - where we have used that the term labelled X in (B11) is non-negative (as (1 − aα) > (1 − a) and γβ < 1) to bound XP 2 ≤ XP , and defined The same inequality can be shown to hold in the other cases, with The same way, we can also obtain for 0 ≤ b ≤ 1 where (or correspondingly -with α and γ exchanged w.r.t. c -when b < 0 or b > 1). Starting from (B1), we thus obtain Conjugating both sides of the equation with ABC, we now immediately see that this is nothing but the condition (B1) for the Hamiltonian h k , with a ij and a ji unaltered, and new feasible points c ij and c ji as given by Eqs. (B13) and (B15). Importantly, in case 0 ≤ a ij ≤ 1 for all i, j, using (B5) this yields a feasible point with for all c ij in (12b), which immediately implies a continuity bound on the gap which only depends on the geometry of the model, but not on the lattice size or the specifics of the model at hand. This can be straightforwardly adapted to the case where some a ij are negative or above 1.
As outlined in the main text, this bound can subsequently used, starting from β = 0, to determine a regime for which the gap can be lower bounded by a uniform ∆ 0 > 0. Clearly, tighter bounds can be obtained by choosing more intermediate values of β.

Appendix C: Symmetrization of the SDP for translational invariant systems
Let us see that for translational-invariant systems, the SDP presented in the main text gives a bound δ which is independent of the system size. For simplicity, we consider 1D systems and the case in which Eq. (15) holds. The most general case, as well as the higher dimensional case, can be generalized from this proof.
Let N denote the system size and let us consider periodic boundary conditions. Consider the SDP (12). Let us assume that we have a feasible point for some choice of the variables (a * ij , c * ij ). Let us show that we can redefine these variables, such as we still have a feasible point with the same value of δ, as follows: It is easy to see that the constraints (12b) are still satisfied. Indeed: For the constraints (12c)-(12d) we get: With the new variablesã ij ,c ij we actually have redundant information in (12b) and it suffices to write down a single condition: It is easy to see now that adding a new particle, i.e., N → N + 1, does not change the SDP (since it will only add redundant information to (C3)). Thus, given that for a TI finite system our method provides a non-trivial bound for the gap, the same bound holds for N → ∞. In this case, we recover the result of the martingale method by Ref. [18]. Note that for higher dimension, instead of a single condition (C3) we would have as many as different types of overlaps within the different Hamiltonian terms.
Appendix D: Bound on the norms of the observables Q λ Let us consider the operators Q λ be defined as in Eq. (21), i.e. Q λ = O † λ P O λ , with O λ as defined in Eq. (18) and P being a Pauli string. We consider the infinite-norm Q λ ∞ and from now on we will drop the subindex for simplicity. Let us write the following inequality for the norm of a product of operators: for A and B being invertible operators, it holds that: Eq. (D1) can be proven by using the submultiplicity of the norm, namely: Here we have used that [κ 1,r , κ 1,s ] = 0, so we can just write the inverse m∈ν(λ) e −βκ1,m −1 = m∈ν(λ) e −βκ1,m −1 .
We can now use that e −βκ1,m −1 = e βκ1,m and that e βκ1,m ∞ ≤ e β . With this, it follows that: Let us now study O † λ : The derivation above implies the following bound: Note that in all the bounds above it is present the cardinality of |ν(λ)|, so let us estimate this value in terms of the size |λ|: Where we have used that |ν j | = O(z 2r1 ), where z is the degree of incidence of the graph at which vertex the particles are placed and r 1 is the radius of the operators κ 1 . With all the derivations above, now we can get a lower bound on Q λ : where we have use that P −1 = P . Now we can use submultiplicity of the norm and Eq. (D4) to get: With this and Eq. (D3) and Eq. (D5) now we get: This means that Q λ decays at most as O(e −2β|λ|z 2r 1 ). If we consider now the terms such that |λ| is constant (or at most goes as O(log N ), with N being the system size) then we get that Q λ decays at most polynomially with the system size.

Appendix E: Analysis of the verification test
In this section, we analyze the verification test for the setup presented in Section IV A, that is, in which the verifier can perform local measurements to the state that the prover has prepared. This case is analogous to the one presented in [28]. We modify it here by imposing further restrictions on the prover: we assume that that he can only perform Pauli measurements, which might be more suitable for experiments carried out with NISQ devices.
We first analyze the probability that, for a given tolerance threshold , the empirical mean h j P computed through the measurements (where the subindex P emphasizes that the verifier measures the state that the prover has prepared) is above such threshold, due to statistical noise produced by finite sampling. This statistical fluctuation depends on the variance of the distribution h j P , which is computed by decomposing the operators in the Pauli basis. This variance is bounded in Appendix E 2.
1. Bounds on the fidelity and the number of measurement rounds Let us denote by h j P the random variable result of estimating the expectation value of h j by measuring the state in the Pauli basis, for a number of samples L j . We define this random variable as h j P = 1 Lj Lj k=1h (k) j , whereh j refers to the random variable associated to h j P computed with a single round (see next section and Eq. (E10) for a formal definition ofh). As before, we assume that L j = L j and that we do not recycle samples, i.e., we always use different samples for different local terms h j . In the same way, we denote by H P the estimator of the energy, which is a random variable whose mean we denote by E [ H P ]. We also define F P = 1 − H P δ , with δ being the lower bound on the spectral gap computed in Section III. By Eq. (26) F P is an estimator of a lower bound on the fidelity between the prover's state and the true ground state. Then, we get the following result: for α ∈ (0, 1) and a threshold parameter > 0, it holds that as long as the number of samples L j satisfies: where N is the number of terms in the Hamiltonian and |λ| is its locality.
Proof. By Eq. (26) we get that: Recall that H P = j h j P . We can assume that the numbers of samples used to estimate any h j P is large enough so the central limit theorem applies. Thus, we can assume that h j P ∼ N (E[ h j P ], Var[ h j P ]) follows a Gaussian distribution, which therefore implies that H P is Gaussian itself. Then, we can bound (E3) by: where the inequality follows from bounding the tails oft he Gaussian distribution and the last equality comes from the fact that Var[ H P ] = j Var[ h j P ] since the random variables h P are independent. Note that Var[ h j P ] = Var[h j ]/L j . We will see in the following section that: By substituting this into the previous equation and imposing the bound (E2) on L j we get: which thus proves Eq. (E1). Note that this equation is similar to Eq. (7) and (15) in [28] and a derivation in terms of Hoeffding bounds as in this paper would also apply here -although the difference is again that we do not measure the eigenenergies of h directly but only the expectation values of strings of Pauli operators.

Bound on the variance of h P
Let us now study the properties of the distribution of h P for a fixed local term h. We first introduce some preliminary definitions, some of which were already introduced in the main text, that will help in the upcoming discussion: • We define a set random variables α j , for j = 1, ..., N . Each variable α j can take the values {x, y, z} with equal probability. We denote the set of this variables by B, B = {α j } N j=1 . Note that each element of B corresponds to the basis in which we choose to measure each qubit independently.
• We denote by J any of the subsets of indices of λ(h) = {j 1 , ..., j |λ| }. Note that the cardinal of J goes from 0, ..., |λ|, where we also include the empty set for convenience. In the same fashion, we denote by B J the subset of random variables α j as defined before with corresponding indices j ∈ J.
• Let s J B J = s j1 αj 1 · · · s j |J| αj |J| be the random variable with possible outputs ±1 resulting from measuring the qubits j ∈ J in the basis α j ∈ B J .
With these preliminaries we can now define the following random variableh: Note that in the definition of the variableh in Eq. (E7) we do not sum over all possible basis choices: this variable actually corresponds to the computation of h P with a single round, in which we only measure the state once in some basis α = (α 1 , ..., α N ). Note that this is different from Eq. (30), in which we averaged over different rounds with the same basis choice (one could in principle also perform an analysis with multi-rounds, but then R(γ) is itself a random variable, which makes the analysis more cumbersome). With this single sample, we can only give an estimator for the following marginals σ j αj , σ j1 αj 1 σ j2 αj 2 , . . . , σ j1 αj 1 · · · σ j |λ| αj |λ| , where α j k ∈ B. Let us now compute the expectation value ofh : E h = o 0 + |λ| j=1 γj =x,y,z 3 · P (α j = γ j )E o j αj s j αj + |λ| j,k=1 γj ,γ k =x,y,z 3 2 · P (α j = γ j , α k = γ k )E o j αj ,α k s j αj s k α k + + . . . + γj 1 ,...,γj |λ| 3 |λ| · P (α j1 = γ j1 , ..., α j |λ| = γ j |λ| )E a Since α j can take the values {x, y, z} with equal probability we have that P (α j = γ j ) = 1 3 , ∀j, ∀γ j . Moreover, the variables α j are independent and therefore the joint probability can be factorized: P (α j1 = γ j1 , ..., α j k = γ j k ) = P (α j1 = γ j1 ) · · · P (α j k = γ j k ) = 1 3 k , for k = 1, ..., |λ|. Note that Eq. (E10) leads to the same expression as Eq.(E8) and therefore: In statistical terms,h is an unbiased estimator for h . We can now compute the variance of the variableh . Now, note that s J B J and s J B J are not independent, since J and J may have non-empty intersection. Therefore we express the variance ofh, defined by the sum in Eq. (E7), in terms of the covariance of its terms: (E14) Note that o 2 max depends on the particular local term that we are considering. We can get rid of this dependence by bounding it by the following value: which follows from the fact that tr h 2 ≤ 2 |λ| .