Bond dimension witnesses and the structure of homogeneous matrix product states

For the past twenty years, Matrix Product States (MPS) have been widely used in solid state physics to approximate the ground state of one-dimensional spin chains. In this paper, we study homogeneous MPS (hMPS), or MPS constructed via site-independent tensors and a boundary condition. Exploiting a connection with the theory of matrix algebras, we derive two structural properties shared by all hMPS, namely: a) there exist local operators which annihilate all hMPS of a given bond dimension; and b) there exist local operators which, when applied over any hMPS of a given bond dimension, decouple (cut) the particles where they act from the spin chain while at the same time join (glue) the two loose ends back again into a hMPS. Armed with these tools, we show how to systematically derive `bond dimension witnesses', or 2-local operators whose expectation value allows us to lower bound the bond dimension of the underlying hMPS. We extend some of these results to the ansatz of Projected Entangled Pairs States (PEPS). As a bonus, we use our insight on the structure of hMPS to: a) derive some theoretical limitations on the use of hMPS and hPEPS for ground state energy computations; b) show how to decrease the complexity and boost the speed of convergence of the semidefinite programming hierarchies described in [Phys. Rev. Lett. 115, 020501 (2015)] for the characterization of finite-dimensional quantum correlations.


I. INTRODUCTION
The study of condensed matter phases depends crucially on our ability to determine the properties of the ground state of local Hamiltonians defined over a lattice. Not only this problem has been shown to be QMA-hard for general Hamiltonians [1,2], but already for mesoscopic systems (of n ∼ 40 particles), even storing the description of a general quantum state in a normal computer becomes an impossible task. This forces condensed matter physicists to resort to quantum state ansatzs in order to understand and study the properties of matter in the low temperature regime.
One ansatz that has proven very useful in this respect is the family of Tensor Network States (TNS) [3], a class of many-body wavefunctions of complexity fixed by a parameter known as bond dimension. In the last few years, TNS have been successfully used to approximate the low energy sector of local Hamiltonians of spin lattices of different dimensions [4][5][6]. The ability to approximately compute expectation values in an efficient manner, together with the possibility to conduct optimizations in the thermodynamical limit [7][8][9][10] makes TNS one of the very few avenues to understand the physics of strongly correlated systems.
Because of all the above and further theoretical considerations, in the last years, almost every talk about TNS starts with the speaker reminding the audience that 'TNS of low bond dimension are the only physical states of condensed matter systems', or, equivalently, that 'all other rays of the Hilbert space of a many body system are not physically realizable'.
Suppose that we take this last claim at face value. That is, we postulate that the laws of Nature are such that the states of condensed matter systems at low temperature are representable via convex combinations TNS of low bond dimension. This is in effect a physical the-ory; as such, its limits must be explored to determine to which degree the theory is falsifiable, and thus scientific. Given that all one can hope to estimate in the lab are the expectation values of certain k-local observables, what makes TNS special, when compared to any other quantum state? How do they transform under local operations? Given experimental data, can we prove that the underlying quantum state cannot possibly have a low bond dimension, i.e., can we falsify a TNS model ?
We lack tools to answer these questions. Note that the naive scheme of lower bounding the bond dimension by estimating the rank of reduced density matrices only works with pure TNS and not convex combinations thereof. Moreover, the physical scenarios which we consider here just allow the experimentalist to estimate averages of two-body reduced density matrices. In addition, the variational methods so commonly used in condensed matter physics to optimize over TNS of fixed bond dimension are useless to refute a TNS model: such a task would require relaxation, rather than variational techniques.
In this paper, we will address these problems for Matrix Product States (MPS) [11], a class of TNS used to model non-critical one-dimensional spin chains.
We start by deriving two surprising features of MPS: first, for any D we identify local operators which annihilate all MPS of bond dimension smaller than or equal to D. Second, for any D we prove the existence of local operators which, when applied over any MPS of bond dimension D, decouple (cut) the particles where they act from the spin chain while at the same time join (glue) the two loose ends back again into a MPS.
Armed with these notions, we will define a family of k-local operators with negative eigenvalues whose expectation values are nonetheless positive for all MPS of a given bond dimension D. Each such operator can be used to certify that the quantum state of a non-critical spin chain does not admit an MPS representation of bond dimension D. Moreover, this partial characterization of the dual cone of MPS allows us to devise general feasibility tests, or automated criteria to falsify MPS models given limited data about the underlying quantum state, such as a number of experimentally available expectation values.
In addition, we will construct instances of local Hamiltonians of arbitrarily many qubits for which a blind application of MPS-based optimization methods would fail to estimate the ground state energy. This construction can be generalized to other TNS for optimizations over spin lattices of higher spatial dimensions. We will also exploit the low dimensionality of the space spanned by MPS and the notion of cut-and-glue operators to decrease the complexity and boost the speed of convergence of the semidefinite programming (SDP) [12] hierarchies described in [13,14] for the characterization of finitedimensional quantum correlations.
The structure of this paper is as follows: first, we will introduce local Hamiltonians and MPS, and also a couple of notions from the theory of Matrix Algebras. Then we will reveal a connection between MPS and polynomials of non-commuting variables, which will allow us to derive non-trivial structural properties of MPS. Next, we will use these properties to explore the limits of MPS models and improve the SDP relaxations proposed in [13,14]. Finally, we will discuss how some of our results generalize to Projected Entangled Pairs States (PEPS) [5].
Before we proceed, though, a disclaimer is in order: long after the completion of this work, we were made aware that the connection between matrix algebras and MPS had already been pointed out by R. Werner in 2006 [15]. In this encyclopaedia article, Werner also observes that the dimensionality of the space spanned by MPS of a fixed bond dimension is polynomial on the system size.

II. MATRIX PRODUCT STATES AND THE THEORY OF MATRIX ALGEBRAS
Consider a quantum system composed of n distinguishable particles, each of which has local dimension d. A general pure state |ψ of this ensemble hence lives in (C d ) ⊗n , and so we require d n complex parameters to describe it.
An n-site Matrix Product State (MPS) is a state of the form where ω, A 1 , ..., A d are D × D matrices. To be precise, in general MPS the matrices {A i } d i=1 are taken to be site-dependent [11]. In the following, though, whenever we refer to MPS, we will have in mind eq. (1). The matrix ω is a boundary condition, while the parameter D is known as the bond dimension of the state. In order to distinguish it from d, the latter is also called the physical dimension.
It can be proven that, for D high enough, all states can be expressed as in eq. (1). However, we will be interested in systems where the value of D does not grow much with the system size n. Note that an MPS of whatever size can be described with O(dD 2 ) complex parameters. Hence, as long as D is not very big, it pays to use this approximation. Finally, notice that, taking ω = I D , the state becomes invariant with respect to the permutation 1 → 2 → ... → n → 1. Such states are called finite translational invariant MPS.
For low values of D, computing expectation values of product operators in an MPS can be carried out in an efficient way. Thus MPS are regularly used to approximate the ground state of k-local Hamiltonians of onedimensional systems, i.e., Hamiltonians of the form: where h j acts non-trivially on the space of the particles j, j + 1, ..., j + k − 1. Given H, we will denote by H D the minimum average value of H achievable with MPS of bond dimension D.
Note that we can always choose 11]. This allows us to perform calculations in the thermodynamic limit, i.e., n → ∞. Indeed, in such a case, the m-site reduced density matrix ρ m of the state under consideration is equal to where σ ≥ 0 satisfies tr(σ) = 1 and i A † i σA i = σ, and the sum runs over all vectors i, j ∈ {1, ..., d} m . With the latter condition, the above is a Translational Invariant (TI) MPS or an infinite MPS (iMPS).
If the state of a finite spin chain can be expressed as a convex combination of MPS of bond dimension D, we will say that it admits an MPS model of bond dimension D. Furthermore, if the TI state of an infinite spin chain can be expressed as a convex combination of iMPS of bond dimension D, we will say that it admits an iMPS model. Unless otherwise specified, whenever we refer to MPS in this paper we will mean states of the form (1) (not necessarily TI).
We now digress momentarily from the topic of MPS to the theory of matrix algebras. A matrix polynomial identity (MPI) F (X) for dimension D is a polynomial of noncommuting variables X 1 , ..., X d that vanishes when evaluated with matrices of dimension D or lower. For example, any commutator [X i , X j ] is a MPI for D = 1. A more elaborate example is [[X 1 , X 2 ] 2 , X 3 ]; this polynomial vanishes when evaluated with 2 × 2 matrices.
MPIs do not only exist for dimensions 1 and 2. In general, it can be proven that any 2D-tuple of D × D matrices X 1 , ..., X 2D must satisfy the standard identity F 2D (X) = 0 [16], where Here S N denotes the set of all permutations π of N elements. It can also be shown that MPIs for matrices of dimension D must necessarily have degree 2D or higher [16].
A concept related to MPIs is that of central matrix polynomials, or polynomials P (X) of noncommuting variables which are proportional to the identity when evaluated with D × D matrices. E.g.: in D = 1, any polynomial can be interpreted as a central polynomial. In D = 2 we already saw an example, namely the polynomial [X 1 , X 2 ] 2 . As with MPIs, it can be proven that non-trivial central polynomials (i.e., central polynomials which are not MPIs) exist for all dimensions D [16].

III. THE PHYSICS OF MPS
We will next establish a relation between matrix polynomials and many-body quantum states. From this link, non-trivial structural properties of MPS will follow almost straightforwardly. Let X ≡ (X 1 , ..., X d ) be any tuple of d noncommuting variables, and let P (X) be a homogeneous polynomial P (X 1 , ..., X d ) = i1,...im p i1,...,im X i1 ...X im of degree m. By |P (X) we will denote the m-particle vector |P (X) = i1,...im p * i1,...,im |i 1 , ..., i m . It is immediate that, for any n-site MPS |ψ of the form (1), applying P (X)| over particles s + 1, ..., s + m leads to That is, we obtain a state similar to an (n − m)-site MPS, but with an 'impurity' in the middle, namely the matrix polynomial P (A). This notion of interacting with physical sites in order to engineer operators at the virtual level is actually the main idea behind measurement-based quantum computing on MPS [17]. Now, let P (X) be a homogeneous MPI for dimension D of degree m. Then, according to (5), P (X)|ψ = 0, i.e., the local operator P (X)| will have the property of annihilating any MPS with bond dimension D or smaller. Conversely, let |P (X) be an m-particle vector with the property of annihilating all MPS of bond dimension D or smaller. Due to our freedom in choosing the boundary condition ω, it is easy to see that P (X) must necessarily be a MPI. We have just established that the local space spanned by MPS is the orthogonal complement of the space of homogeneous MPIs for dimension D.
The next question to answer is how big these two : Calling a the entries of the matrices A, we thus have that where the sum runs over all monomials w of degree m in a and degree 1 on the entries of ω, and the vectors {|φ w } w do not depend on the particular values of A, ω (e.g.: the vector corresponding to the monomial ω 11 (A 1 ) m 11 is |1 ⊗m ). The above decomposition allows us to bound the dimensionality of H M P S D,m simply by counting the number of such monomials. The result is where the D 2 factor stems from the number of entries of the boundary condition ω. Relation (7) implies that the dimension of H M P S D,m increases polynomially with the system size m, contrarily to the total local space dimension, which increases as d m . In the limit of high m, the space of MPIs is therefore exponentially bigger than H M P S D,m . Both the identification of H M P S D,m with the orthogonal complement of H M P I D,m and the polynomial bound on dim(H M P S D,m ) appear in Werner's encyclopaedic article [15].
In Appendix A we describe two efficient algorithmic procedures to generate an orthonormal basis for H M P S D,m . This allows us to ascertain the exact dimensionality of the spaces H M P S D,m , for whatever values of D, d, m. Some results are presented in Table I.
Let h ∈ B(H M P I D,m ) be a self-adjoint operator acting on the space of MPIs of dimension D, and suppose that we integrate it in an n-body k-local Hamiltonian H, with k ≥ m. That is, suppose that the Hamiltonian of the system is H = H + h, with H given by eq. (2).
The above discussion implies that MPS of bond dimension D or lower 'will not see' such a term, i.e., H |ψ = H|ψ for all MPS |ψ = |ψ(ω, A, n) , with ω, A ⊂ B(C D ), n ≥ m. Elaborating on this, we find The result also holds when we restrict the Hamiltonian minimization to finite TI MPS (with ω = I D ). Moreover, it can be extended to iMPS and TI Hamiltonians in the infinite spin chain. The Hamiltonian can also be taken O(D)-local at the cost of increasing the physical dimension of the particles.

See Appendix B for a proof.
A blind application of the proverbial method of minimizing a Hamiltonian via MPS of increasing bond dimension until the sequence of energy values appears to converge hence risks getting stuck at a suboptimal point. Admittedly, the limitations implied by the Proposition do not pose a practical threat for usual studies of onedimensional non-critical chains, since current MPS-based algorithms allow reaching bond dimensions of order 100 in a normal computer (way beyond the locality of Hamiltonians of physical interest). This is no longer the case, though, for some other classes of TNS, for which we lack good optimization schemes and Proposition 1 also extends, see below.
Suppose now that we choose P (X) in eq. (5) to be a central polynomial for dimension D. Then P (A) = p(A)I, where p(A) is a scalar. The state |ψ will hence get projected into the state p(A)|ψ , with This is again a MPS with the same boundary condition ω and matrices A 1 , ..., A d , but where m particles are just missing. Now, divide the space of m-degree homogeneous central polynomials into classes [P ] = P + M P I, i.e., two central polynomials P 1 , P 2 belong to the same class if and only if P 1 − P 2 is an MPI. These classes form a vector space Q D,m , the quotient space of m-degree homogeneous central polynomials by MPIs. Let {P i (X)} i be a basis for Q D,m and define the m-local operator C ≡ i |ϕ i P i (X)|, where {|ϕ i } i is any orthonormal set of m-particle states. The effect of C over any MPS of bond dimension D or smaller is to project the m particles where it acts into the pure state |ϕ ≡ i p i (A)|ϕ i , while the remaining particles end up in the state (9). It hence 'cuts' particles s + 1, ..., s + m off the chain and 'glues' the two ends back, see Figure 1.
This sort of operators will be called cut-and-glue operators. Note that, after tracing out the particles cut, a cut-and-glue operator acting over a TI MPS leaves the spin chain in the same quantum state (modulo normalization).
In Appendix C we sketch efficient procedures to generate a basis for Q D,m . Table II gives an idea of how the dimensionality of Q D,m scales with the system size and the bond dimension in qubit ensembles. Surprisingly, it turns out that the dimensionality of Q 2,m in Table II follows the sequence of coefficients in the power series expansion of the Poincaré series P (C 2,2 ; t), which is sequence A096338 in the On-Line Encyclopedia of Integer Sequences [18]. We conjecture that the above observation holds true for entries beyond m = 15 in Table II as well.
Both annihilation and cut-and-glue operators are important structural features of MPS. Unfortunately, even for D = 2 their implementation in the lab would require the ability to switch on non-trivial four-interaction terms. This is experimentally challenging, given that in many experimental setups only 2-local operators are accessible. Fortunately, there is a cleverer way to exploit our findings.
The notions of anihilation and cut-and-glue operators allow us to define a family of local operators h whose average value is non-negative when computed with nsite MPS of bond dimension D or smaller. Call P the projector onto the space H M P S D,n , and consider all m-local operators h which satisfy: where f ≥ 0, C j is a cut-and-glue operator acting nontrivially over particles 1, ..., j and g j is an entanglement witness [19] with respect to the partition 1, ..., j|j+1, ..., n (namely, g j ≥ 0 for all quantum states separable with respect to the said partition). Clearly, h D ≥ 0.
Given an arbitrary (k-local) operator H, consider the problem of maximizing µ ∈ R such that H − µ admits a decomposition of the form (10). Then, for any feasible µ, H D ≥ µ. If the minimum eigenvalue of H happens to be smaller than µ, then H − µ can be regarded as a bond dimension witness: an expectation value for H below µ would prove that the underlying quantum state of the system does not admit a MPS model of bond dimension D.
Regretfully, the optimization proposed above requires an implicit characterization of entanglement witnesses, a problem known to be NP-hard [20]. A way out is to simply demand g j to belong to a class of entanglement witnesses which are easy to describe. An obvious choice is the set of all operators which are Positive (semidefinite) under Partial Transposition (PPT) [21]. With this restriction on g j , the maximization of µ can be formulated as a semidefinite program (SDP), a class of convex optimization problems which can be solved efficiently [12].
The dual of this program would be an optimization over all quantum states ρ ∈ B(H M P S D,m ) such that, for all j, C j ρC † j is PPT for the partition 1, ..., j|j + 1, ..., n. Since H M P S D,m grows polynomially with m, for small D a normal computer can reach large values of m. Moreover, playing with the displacement operator, it is easy to derive a hierarchy of SDPs for the characterization of iMPS models.
Rather than describing these tools in detail [the reader can find a full description of the SDP programs in Appendix D], we will illustrate how these methods work with a practical example. Consider an N -site spin 1/2 chain, and suppose that, via neutron interferometry, we estimate the expectation value of the XXX Heisenberg Hamiltonian We wonder whether our experiment can be explained with an iMPS model of low bond dimension. Take the number of particles in the chain to be small, say N = 7. We are interested in determining whether our system admits an MPS model of low bond dimension. For low N , the minimum eigenvalue of H N can be computed exactly, and so we find that the minimum average energy per interaction term is E ≡ min 1 6 H 7 ≈ −0.4727. On the other hand, an SDP optimization over 7-site normalized density matrices ρ ∈ B(H M P S 2,7 ) satisfying C j ρC † j , PPT for j = 5, 6 returns the greater value E 2 ≡ 1 6 H 7 2 ≥ −0.4065. This optimization, and all subsequent ones, was carried out with the SDP solver MOSEK [22]. The XXX Heisenberg Hamiltonian can thus be interpreted as a displaced bond dimension witness: an expectation value smaller than -0.4065, within reach given the lower value of E, would signify that the state of the spin chain cannot have a MPS model of bond dimension D = 2.
Refuting D = 3 MPS models for N = 7 is impossible with the tools developed so far, since H M P S  Tables I and II. To make matters worse, SDP optimizations for N ≥ 9 are too memory-demanding for a normal computer. However, if we drop the PPT condition, the resulting SDP can be seen equivalent to projecting 1 N −1 H N on the subspace H M P S 3,N and finding the minimum eigenvalue of the resulting operator. This simplified method allows us to reach greater values of N , at the price of losing robustness in our bounds. With this trick, for N = 13 we obtain a bound E 3 ≡ 1 12 H 13 3 ≥ −0.44958, slightly bigger than the minimum energy density E = −0.46044 achievable.
For N 1, we can take the system to be approximately translational invariant, so this time we want to refute iMPS models of low bond dimension for our system. For D = d = 2 iMPS models, it can be shown that spin chains are symmetric under parity, and hence satisfy non-trivial linear constraints. However, if our experimental setup does not allow us to estimate quantities of the sort , with t = s, we must again rely on inequalities, in which case the XXX Heisenberg Hamiltonian can also serve as a witness. It is a standard result (see, e.g., [23]) that the ground state energy density E = lim N →∞ Now, how to derive bounds for iMPS models? A possibility would be to compute the minimum expectation value of 1 N −1 H N for high N via the MPS SDP relaxation used above. Intuitively, increasing values of N should give better and better approximations to the optimal iMPS value for the energy density. We chose, though, to use the slightly better approximation of optimizing the value of 1 4 σ 1 · σ 2 over reduced density matrices subject to the constraints above and the extra condition ρ ∈ S D,N , where S D,N denotes the span of the N -site reduced density matrices of iMPS of bond dimension D. This space can be characterized using similar techniques as the ones we applied to compute H M P S D,m . Optimizing over 8-site normalized reduced density matrices ρ ∈ B(H M P S 2,8 ) satisfying C j ρC † j , PPT for j = 5, 6, 7 and ρ ∈ S 2,8 we find that E 2 ≡ lim N →∞ 1 N −1 H N 2 ≥ −0.3378. An average energy lower than the last value will hence refute all iMPS models of bond dimension D = 2.
This example is very illuminating in that it allows to appreciate the relevance of cut-and-glue operators for this class of optimizations. For, if we drop the PPT conditions above, the lower bound on E 2 output by the computer decreases to −0.4246. This is still bigger than E, and so it also defines a bond dimension witness. However, it is one order of magnitude less robust than the previous one, and so its violation is more challenging from an experimental point of view.
As a final example, we consider the Majumdar-Ghosh Hamiltonian [24]: The expectation value of this operator can also be estimated experimentally via neutron difusion. In the thermodynamical limit N → ∞, the minimum energy density of H M G is E = − 3 8 = −0.375, achievable with an iMPS of bond dimension D = 3 [11]. In contrast, an 8 th -order SDP relaxation over iMPS of bond dimension D = 2 gives E 2 ≥ −0.2593. We have just derived a bond dimension witness with a large gap between iMPS models with D = 2 and D = 3. On the negative side, though, our lower bound for E 2 is significantly lower than the best upper bound E 2 ≤ −0.125 we found using the Amoeba variational method [25].
Note that the former SDP methods can be easily turned into feasibility tests. Indeed, determining the existence of a state with the properties above compatible with some partial information we may hold about the quantum state of the spin chain (such as, e.g., the average value of a number of 2-local observables) can also be cast as an SDP. This procedure can help an experimentalist to refute the existence of an MPS or iMPS model for the state he/she prepared in the lab, without the need of guessing the 'right' bond dimension witness to do the job.
An immediate question is whether the SDP hierarchy of relaxations for iMPS models sketched above is complete, in the sense that it allows us to detect any state lacking an iMPS model by taking N sufficiently large (the SDP relaxation for general MPS is clearly not complete). In this regard, notice that H M P S 1,N corresponds to the symmetric space of N particles, and Q 1,N = H M P S 1,N . Hence for D = 1 the hierarchy reduces to just imposing that the overall state is symmetric and PPT with respect to any bipartition. This is actually the Doherty-Parrilo-Spedalieri (DPS) method for entanglement detection [26], and convergence follows from the quantum de Finetti theorem [27]. Similarly, H M P S ∞,N = C d N , Q ∞,N = {0} and S ∞,N is the span of all TI states. For D = ∞ the hierarchy is therefore computing the minimum expectation value of a Hamiltonian term over N -site states whose N −1-site reduced density matrices coincide whenever we remove the first or the last site. This is essentially a reformulation of Anderson's approximation [28], where convergence is also known to hold. One would be tempted to claim that our SDP hierarchy should converge as well for all intermediate values of D, but we leave this matter open.

IV. APPLICATIONS FOR OPTIMIZATIONS OVER FINITE DIMENSIONAL QUANTUM CORRELATIONS
In [13,14], a hierarchy of SDP relaxations is presented to characterize the statistics of finite-dimensional quantum systems. This hierarchy relies on the notion of moment matrices. Given a quantum system in state σ ∈ B(C D ), with (self-adjoint) operators X 1 , ..., X d−1 ⊂ B(C D ), its n th order moment matrix M is a matrix whose rows and columns are labeled by monomials u of X 1 , ..., X d−1 of degree smaller than or equal to n, with entries given by M u,v = tr(u † (X)σv(X)). In [13,14], it is proposed to relax the requirement of M admitting a quantum representation by demanding M ≥ 0 and M ∈ M D , where M D denotes the space spanned by moment matrices with quantum representations of dimension D.
A disadvantage of this method is that, for fixed d, D, the complexity of implementing the hierarchy increases exponentially with the index n of the relaxation. In the following, we show that every n th -order moment matrix with a quantum representation of dimension D can be interpreted as a conic combination of n-site MPS with bond dimension D. This will allow us to devise an MPS-based algorithm that carries out an improved version of the n thorder relaxation described in [13,14] in time polynomial in n.
Define X d ≡ I D , and consider vectors of d indices with values in {1, ..., d}. Then, for any index vector i ∈ {1, ..., d} k , we can associate the monomial u(X) i ≡ X i1 ...X in . This procedure gives an over-representation of the set of monomials of degree smaller than or equal to n. Now suppose that, by repeating rows and columns, we enlarge the n th -order moment matrix M of the system to an d k × d k matrixM such that i|M | j = M u i ,v j , with | i = |i 1 ...|i n and similarly for | j . The 'enhanced' moment matrixM can then be written as This is just the transpose of a conic combination of MPS of bond dimension D; much like eq. (3), but without the condition i A † i σA i = σ. As such, its support is contained in H M P S D,n ; more precisely, in the analog set for MPS with X d = I plus any other extra restriction in the variables X. HenceM can be fully specified by a number of parameters polynomial in n.
Most interestingly, one can compute special cut-andglue operators C for this kind of MPS. The convergence of the scheme can therefore be boosted by demanding extra positive semidefinite constraints such as (C jM C † j ) Tj ≥ 0. 'Localizing matrices' of the form M q u,v = i, j tr(X † i σX j q(X))| i j|, defined in [14] to model semi-algebraic conditions of the sort q(X) ≥ 0, can be treated in a similar way.

V. EXTENSION TO PROJECTED ENTANGLED PAIRS STATES
MPS can be understood as elements of a larger class of TNS called Projected Entangled Pairs States (PEPS) [5]. Such states are used to approximate the low energy sector of local Hamiltonians describing particles sited in square lattices of arbitrary spatial dimensions. While MPS are defined via tensors with one physical index (i = 1, ..., d) and two bond indices (the column and row indices of the matrices A 1 , ..., A d ), N -dimensional PEPS are defined via contractions of tensors with 1 physical index and 2N bond indices. MPS can therefore be regarded as one-dimensional PEPS. For illustration, in Figure 2 the tensors of a two-dimensional PEPS are represented by circles, while physical [bond] indices are denoted by red [orange] lines. It is natural to ask whether some of the structural features we derived for MPS also extend to PEPS of higher spatial dimension.
Consider the space spanned by PEPS of bond dimension D and physical dimension d in a given region R of the lattice, with boundary ∂R, see Figure 2. As with MPS, we can express any PEPS in R as in (6), where each monomial w has degree 1 on the boundary condition and degree |R| on the tensor A generating the PEPS. The local space spanned by PEPS is hence bounded by D |∂R| poly(|R|). This bound must be compared with the total dimensionality of the physical space in R, namely, d |R| . Provided that d |R| > poly(|R|)D |∂R| , we will find non-trivial operators in R which will annihilate all PEPS with bond dimension D or smaller. If the spatial dimension of the lattice is N , taking R to be a hypercube of size L, with volume |R| = L N and surface area |∂R| = 2N L N −1 , this will happen for high enough L.
We have just proven the existence of tensor polynomial identities, i.e., local vectors |φ which annihilate all PEPS of bond dimension D. Given that generic PEPS are the ground state of a unique parent Hamiltonian, it is easy to prove a weaker version of Proposition 1 for PEPS of arbitrary spatial dimension. Namely, that the chain of identities will break at some point beyond H D (not necessarily at D + 1). Let us remark that, contrarily to MPS variational algorithms, current tools for In such a predicament we may be eager to believe that the last estimation is a good approximation to the ground state energy, if the corresponding optimizations over lower bond dimensions returned similar results... and we could be wrong, as the arguments above show.

VI. CONCLUSION
In this work, we have presented two highly non-trivial structural properties of MPS, namely, the existence of annihilation and cut-and-glue operators. We used these notions to prove several results concerning the limitations of the MPS approximation. Along the way, we raised a number of important open questions.
First, it would be desirable to find closed formulas for the dimensionalities of H M P S D,m and Q D,m . Even though we have efficient methods to calculate these exactly, large values of m require a considerable amount of computational time (hence the missing entries in Tables I, II). Perhaps the connection with the Poincaré series P (C 2,2 ; t) can be exploited in this regard.
Another interesting problem is whether our SDP hierarchy of relaxations to refute iMPS models is complete or can be further improved.
Regarding completeness, for D = 1 the proof of convergence follows from the quantum de Finetti theorem [27]. For D = ∞, the hierarchy is just a reformulation of Anderson's approximation [28], whose convergence was established long ago. A convergence proof for all other values of D would not only provide us with an alternative definition of iMPS, but most likely would involve an intermediate result of depth comparable to the quantum de Finetti theorem.
As for improvement, a promising avenue to boost the speed of convergence of the hierarchy is to incorporate to our codes entropic constraints of the form S(ρ 1,...,k ) ≥ S(ρ 1,...,k−1 ) for 2 ≤ k ≤ n, as in [29]. These hold for any TI state; and, in particular, for iMPS. Although not reducible to SDPs, the corresponding problems can nonetheless be attacked with the tools of convex optimization theory. Considering the space spanned by several copies of a MPS should also help.
Admittedly, it is difficult to believe that these ideas can ever lead to non-trivial restrictions for MPS models with D > 4. Devising new tools for the characterization of MPS models of high bond dimension is hence an important matter.
Finally, it is intriguing whether the analogs of cut-andglue operators for MPS also exist for PEPS of higher dimension. The action of such local operators over an arbitrary PEPS would be to project the region where they act into a pure state and interconnect the links of the particles in the boundary. Appropriately tamed, such operators would allow transforming tensor network states of different type into each other by means of fixed (i.e., state-independent) local operations, very much like graph states transform into each other [30]. In this case, however, computational explorations face the exponential complexity of characterizing the space spanned by PEPS of dimensions two and higher.
Determining the exact dimensionality of H M P S D,m and deriving an orthonormal basis for this subspace can be done via two different procedures. First, given an arbitrary weight w(ω, A) > 0, we will call polynomial MPS any m-site states of the form where f (ω, A) is a homogeneous polynomial of the components of ω (with degree 1) and A (with degree m). A comfortable possibility is to take ω, A real and the trivial weight w(ω, A) = 1. Defining t ≡ dim(H M P S D,m ), our task is to find a set of polynomials {f i } t i=1 such that f M P S i |f M P S j = δ ij and span{|f M P S i } = H M P S D,m . This can be seen equivalent to diagonalizing the kernel and taking the eigenvectors (polynomials) with non-zero eigenvalue.
Alternatively, we can simply sequentially generate real random D×D matrices ω j , A j 1 , ..., A j d and use them to define the sequence of random MPS (|ψ(ω j , A j , m) ) j . Exploiting the fact that the overlap between two MPS can be computed efficiently, one can apply a Gram-Schmidt process to the previous sequence of MPS, thus obtaining an orthonormal basis for H M P S D,m , whose elements are finite linear combinations of MPS. The goal of this section is to prove Proposition 1 in the main text, which reads: For the proof we need two intermediate results, namely: 1. There exist TI O(D)-local Hamiltonians of arbitrarily many particles whose unique ground state is a MPS of bond dimension D. This is proven in Section B 1.
2. For any D > 1, there exists a bivariate noncommutative homogeneous polynomial F (X 1 , X 2 ) of degree O(D 2 ) that is a MPI for matrices of size D − 1 × D − 1, but not for matrices of size D × D. This will be proven in Section B 2.
These two results will be combined to demonstrate the Proposition in Section B 3.

Parent Hamiltonians and MPIs
A set of matrices A 1 , ..., A d ∈ B(C D ) satisfies the injectivity condition if there exists k such that the products {A i1 ...A i k } span all of B(C D ). From [11], we know that any (TI) n-site MPS whose matrices satisfy the injectivity condition for some order k can be seen as the unique ground state of a (TI) 2k-local Hamiltonian (provided that n ≥ 2k).
We will next prove that, for any D, there exist matrices B 1 , B 2 ∈ B(C D ) which satisfy the injectivity condition for k = O(D). With the above, this will imply that, for any bond and physical dimensions D, d and any system size n, there exists an n-site TI MPS which arises as the unique ground state of a O(D)-local TI Hamiltonian.
Let d = 2, and consider the matrices |i j| (B2) Note that we can express the projectors {|i i|} D i=1 as linear combinations of {B p 1 : p = 1, ..., D}. Since B s 2 = B 2 for any s ≥ 1, this implies that linear combinations of the (degree 2D + 1) products B p 1 B 2D−p−q+1 2 B q 1 can generate the matrices |i i|B 2 |j j| = 1 D |i j|, which span B(C D ).

MPIs for dimension D − 1 which cease to be identities in dimension D
In this section we will prove that, for any D, there exists a polynomial F (X 1 , X 2 ) of degree O(D 2 ) which is a MPI for dimension D − 1, but not for dimension D.
Note that the matrices {P j (B)} j have the peculiarity that the only permutation of them which does not vanish is P 1 (B)P 2 (B)...P 2(D−1) (B).

Putting all together
Now we are ready to prove Proposition 1.
Proof. Take a TI n-qubit MPS |ψ with bond dimension D , and build its TI parent Hamiltonian H [11]. Such is a k-local operator with the properties H ≥ 0, H|ψ = 0 and |ψ being the only ground state of H. From Section B 1 we know that |ψ can be chosen injective and such that H has interaction strength k = O(D ). Since |ψ is injective, it cannot be expressed as an MPS of bond dimension D − 1 (because, e.g., |ψ has Schmidt rank greater than D − 1). It follows that H Note that, if we are entitled to play with the physical dimension of the system, we do not need to invoke Section B 2 at all. Indeed, it suffices to set d = 2D and take h to be the standard identity for dimension D (which, having degree 2D, cannot be an MPI for dimension D + 1). The corresponding family of Hamiltonians H λ would then be O(D)-local, rather than O(D 2 )-local. polynomials are central in D = 1). The method hence reduces to the Doherty-Parrilo-Spedalieri (DPS) method for entanglement detection [26]. The DPS method approximates the set of states of the form dφp(φ)|φ φ| ⊗m , with p(φ) ≥ 0 by partial traces of the set of n-symmetric states positive under partial transposition. It can be shown to converge by virtue of the quantum de-Finetti theorem [27]. For D = ∞, program (D2) is equivalent to computing the (n − m) th Anderson bound for the Hamiltonian H [28], and convergence can be proven easily. It is an open question whether the hierarchy (D2) converges for other values of D.