Certifying ground-state properties of many-body systems

current implementation, neither variational nor relaxation methods offer provable bound on other observables in the ground state beyond the energy. In this work, we show that the combination of the two classes of approaches can be used to derive certifiable bounds on the value of any observable in the ground state, such as correlation functions of arbitrary order, structure factors, or order parameters. We illustrate the power of this approach in paradigmatic examples of 1D and 2D spin-one-half Heisenberg models. To improve the scalability of the method, we exploit the symmetries and sparsity of the considered systems to reach sizes of hundreds of particles at much higher precision than previous works. Our analysis therefore shows how to obtain certifiable bounds on many-body ground-state properties beyond energy in a scalable way.


I. INTRODUCTION
The quantitative description of many-body quantum systems is one of the most important challenges in physics.A standard formulation of the problem consists of N particles each described by a Hilbert space of dimension d with interactions encapsulated by a Hamiltonian H acting on C d N .Typical questions of interest are the study of the system Hamiltonian evolution, the computation of the energy spectrum, or the characterisation of thermal and ground-state properties.A brute force approach to these problems requires diagonalisation of the Hamiltonian, or more generally, dealing with matrices whose dimension (equal to d N ) grows exponentially with the number of particles.It is therefore intractable beyond small clusters of particles.
Among quantum many-body problems, the study of ground states plays a central role due to its relevance for the understanding of the low-energy phases in the system, and in particular for the study of genuine quantum correlation properties without a classical analog [1].Formally, a ground state |ψ GS ⟩ of a quantum Hamiltonian H is a state of minimal energy, that is, a minimiser to the following problem: As mentioned, if |ψ⟩ is decomposed in some given basis of the Hilbert space with an exponentially-growing number of parameters, the exact solution to this optimisation problem quickly becomes intractable when increasing the system size.In fact, the very enumeration of the exact ground-state coefficients in a given basis is out of reach.
To solve this issue, the standard approach consists of finding approximations to the optimisation that provide a much better scaling in terms of computation.By far, the most popular approach is given by variational methods.There, the minimisation in Eq. ( 1) is restricted to a subset of Ansatz states A for which the computation of the expectation value ⟨ψ|H|ψ⟩ and its minimisation is scalable with the number of particles, From an optimal solution state |ψ A ⟩, one then computes the value of some physically-relevant observables o A = ⟨ψ A |O|ψ A ⟩.In these methods, the hope is that the set of Ansatz states A is suitably chosen so that the obtained energy and state are close to the unknown exact values, E A ∼ E GS and |ψ A ⟩ ∼ |ψ GS ⟩, and so are other physically-relevant quantities, o A ∼ o GS = ⟨ψ GS |O|ψ GS ⟩.Variational methods suffer two major limitations: First, they provide no guarantee that the derived upper bound E A is close to the exact ground-state energy E GS .Second, and more problematically, even with a promise that E A is close to E GS , there is absolutely no guarantee that the state |ψ A ⟩ is close to the ground state |ψ GS ⟩ (unless one has more information about the system of interest, such as its energy gap).Hence, for observables other than the energy, it is not known how the computed value o A relates to the actual ground-state value o GS , and in particular whether o A represents a lower or upper bound to o GS .As it turns out, o A may significantly differ from o GS , as for instance strikingly observed in some fermionic Hubbard models [2].These issues shall be discussed in more details in the core of the paper.
Complementary approaches to variational methods, and which are the focus of the present work, are socalled relaxations of the ground-state problem.The general idea is to minimise the energy over a set of parameters that contains all the physically-possible expectation values ⟨ψ|H|ψ⟩, but also other values that are not allowed by quantum mechanics.As for variational methods, the optimisation in the relaxation has a much better scaling than exact diagonalisation and can be computed for larger system sizes.An example of this approach is given by the semidefinite programming (SDP) relaxation to non-commutative polynomial optimisation problems, which have been considered in many different groundstate problems, see for instance [3][4][5][6], and formalised in [7,8], see also [9,10].This relaxation, also known as the Navascués-Pironio-Acin (NPA) hierarchy, plays a fundamental role in this article and is detailed below.By construction, as the minimisation is performed over a set of solutions that contains all the physical values (and possibly extra values), any relaxation provides a lower bound of the exact ground-state energy E R ≤ E GS .Relaxations can also provide estimates o R to the expectation value of observables O in the ground state, but in contrast to variational methods, one cannot guarantee that these estimates are compatible with some underlying quantum state.Furthermore, as for variational methods, there is no a priori guarantee that these values are close to those of the ground state, o R ∼ o GS , neither whether they represent upper or lower bounds to the ground-state values.In summary, combining the present techniques, all what can be certified about ground states is that its energy lies within the range In this work, we show that the combination of variational methods together with the NPA hierarchy is much richer than previously envisioned, and allows for deriving certified upper and lower bounds on the values of arbitrary observables in the ground state, o GS ∈ [o LB , o UB ].We achieve this by assisting the non-commutative polynomial relaxations with some available upper bound of the ground-state energy as given by variational methods, an approach also considered in [11].This allows for computing lower and upper bounds to any observable that can be expressed as a polynomial in a family of basic observables.Examples of these operators are correlation functions of arbitrary order, and structure factors which characterize long-range fluctuations in many-body systems.We apply this approach to several paradigmatic spin models.We focus on Heisenberg models with local interactions and translation symmetry in one and two spatial dimensions, and exploit these properties to construct SDP relaxations for systems of up to a hundred of particles, obtaining much better bounds to ground-state observables than achieved in previous works.Our approach therefore provides a scalable way to derive provable bounds on ground-state properties, beyond the energy.
The structure of the article is as follows.In Section II we introduce the quantum ground-state problem.We review both variational methods and relaxations, with an emphasis on SDP relaxations of non-commutative polynomial optimisation because of their central role in our analysis, and see how they provide, respectively, upper and lower bounds to the ground-state energy.In Section III we present our main idea and show how to derive certifiable bounds on the ground-state value of any polynomial observable when combining the two approaches.In Section IV we provide several applications and illustrations of our construction.We first introduce the general form of the considered models, and discuss how to exploit their symmetries and the sparsity of the Hamiltonian to reach systems of hundreds of particles.We then apply the method to different Heisenberg models in one-and two-dimensional lattices.In Section V we finally discuss our findings and display our conclusions.

II. GROUND-STATE PROBLEM
We consider quantum systems composed of N particles, whose interactions are described by a Hamiltonian operator.In what follows, and for simplicity, we are going to focus our discussion on systems of finite dimension d on a lattice, although the techniques we discuss also apply beyond this scenario, e.g. to boson or fermion models.We consider systems with local interactions and translational symmetries so that it is possible to define a Hamiltonian for an arbitrary number of parties H N given by the sum of different tensor product terms h i acting on a subset of neighbouring particles, H N = i h i .One is then interested in determining the value o Exactly solving the ground-state problem is computationally too costly already for systems of several tens of particles, as the dimension of the systems grows exponentially with the number of particles, d N .Hence one should abandon looking for an exact solution to the problem and adopt approximations to it, such as variational methods or relaxations, which offer a much more favorable scaling with N .

A. Variational methods
Variational methods restrict the ground-state optimisation to a subset of Ansatz states A. It is then demanded that the number of parameters needed to specify an Ansatz state, |ψ A ⟩, scales polynomially with the number of particles, N , and that the computation of mean values of the operators in the Hamiltonian, ⟨ψ A |h i |ψ A ⟩, is efficient.This allows for solving the energy minimisation over Ansatz states, Eq. ( 2), for systems much larger than those for which an exact diagonalisation is possible.
Mean field is one of the simplest instances of a variational method, where the set of Ansatz states is defined by product states.Here, the number of parameters scales linearly with the system size, N d, and the mean value of the local terms h i is easy to compute.The Density-Matrix-Renormalisation-Group (DMRG) approach has represented a breakthrough in the design of variational methods, for it often allows one to obtain good approximations to the ground-state energy of gapped 1D systems [12,13].It is now well understood that the Ansatz states relevant in DMRG are the so-called Matrix-Product States (MPS), whose description requires N χ 2 parameters [14][15][16].Here χ is the so-called bond dimension that determines the entanglement properties of the MPS state.In fact, product states, used in mean-field calculations, are MPS of bond dimension χ = 1.It is known that DMRG works well for 1D systems because ground states of gapped 1D systems with local interactions can be approximated by MPS of fixed bond dimension, that is, DMRG is optimising over a set of states that contains a very good approximation to the unknown ground state [17].Using insights from entanglement theory, it was possible to generalise MPS to other subset of states, such as Projected-Entangled-Pair State (PEPS) [18] or Multi-Entanglement-Renormalization-Ansatz (MERA) [19][20][21], which may be viewed as special instances of the more general set of Tensor-Network states [17].
Other popular variational states which are not based on tensor networks include resonating valence bond states, introduced in the context of quantum magnetism [22], neural-network quantum states [23], or correlatedplaquette states [24][25][26].
It is not our purpose to review here all variational methods; instead, we emphasize that despite all their remarkable applications in the study of many-body quan-tum systems, variational methods do have intrinsic limitations.First of all, for many systems, it is not known whether the actual ground state can be well approximated by a state in the chosen set of Ansatz states.For instance, mean-field energy values can easily be computed for very large sizes, but it is expected that ground states of generic Hamiltonians are in fact entangled, so that all entanglement properties -which sometimes form defining properties of the phase as in topological quantum matter [27,28] -are inaccessible by construction.Second, even when Ansatz states properly approximate the ground state, the minimisation of the energy may remain hard.While an efficient algorithm exists for ground states of 1D gapped local Hamiltonians [29], it is expected that this is an exception.For instance, ground states of 2D gapped systems are known to be well approximated by PEPS, but the computation of expectation values of product operators with PEPS, including even the norm of the state, is #P-hard in the number of tensors defining the state [30].Third, even if the ground state can be approximated by a given Ansatz state and the computation of the energy is scalable, its minimisation often presents many local minima and therefore one can never guarantee that a good approximation to the ground state has been achieved, E A ∼ E GS .But most importantly, the situation is even worse for other relevant quantities computed from the Ansatz state resulting from the minimisation, as it is completely unknown how they compare to the actual values in the ground state.In fact, there exist simple paradigmatic models displaying a very complex low-energy landscape, with states very close in energy to the ground state, but with significantly different predictions for other quantities -the fermionic Hubbard model being a prominent example [2].
In summary, variational methods have proven extremely useful to analyse ground-state problems; yet strictly speaking all what can be certified is an upper bound of the ground-state energy, E A ≥ E GS .

B. Relaxations
Relaxations of the ground-state problem represent a complementary approach to variational methods.Without again attempting to be fully general, the main idea is rather simple: for a Hamiltonian H N = i h i , the ground-state problem, Eq. ( 1), can be reformulated as: where A simple way to compute lower bounds to ground-state energies of a large system of size N consists of performing exact diagonalisation of the same model for the largest possible size, say K < N .This approach is known as the Anderson's bound [32] and easily follows from the fact that the set of moments compatible with an N -particle quantum state is strictly included in the one compatible with a quantum state of K < N particles (see Appendix B for further details).
Relaxations of polynomial optimisation problems based on SDP provide another way to derive lower bounds.Because of the central role they play in our analysis, we explain the main idea in what follows, while more details about their implementation are given below.These relaxations apply to any optimisation problem of the form [8]: such that: where p, g j and h k are polynomials defined over a set of n operators X = (X 1 , . . ., X n ), the minimisation runs over all possible Hilbert spaces, and where ⪰ refers to operator positivity.For simplicity in the notation, we restrict the explanation to self-adjoint operators, although all what follows also applies to general operators.The solution to this problem may be hard.For instance, deciding whether the solution to this optimisation problem satisfies either p min < 0 or p min > 1 is Turing undecidable [33].
In Refs.[7,8], it was shown how to construct an infinite hierarchy, known as NPA, of monotonically increasing lower bounds to the solution of the problem, p (1) ≤ p (2) ≤ • • • ≤ p (∞) ≤ p min .Under some mild assumptions on the operators, that include the situation in which all the operators X i are bounded, the NPA hierarchy is convergent, that is p (∞) = p min .The main advantage of the hierarchy is that the computation of each lower bound defines a SDP instance and therefore can efficiently be performed.However, the size of the operators involved in each SDP step of the hierarchy is also monotonically growing and, in the limit, involves operators of infinite size.Interestingly, for some problems, convergence is attained at a finite step in the hierarchy, while for others, first steps of the NPA hierarchy provide a good enough bound of the actual solution.
Because of the generality of the formalism, many ground-state problems can be phrased in the language of polynomial optimisation.For instance, fermionic Hamiltonians are often given by polynomials of creation and annihilation operators, which satisfy the anti-commutation relations in form of other polynomials.For spin-onehalf systems, Hamiltonians are defined by polynomials of Pauli matrices on each site.Pauli matrices can be characterised by their algebra while operators on different sites commute, all these constraints having the form of polynomials (note that for the Pauli algebra, these constraints impose that the solution Hilbert space of Eq. ( 4) is of finite dimension).The method therefore provides a rather versatile approach to derive an asymptotically convergent sequence of lower bounds to the ground-state energy of many models of interest, In fact, as already mentioned, it has already been applied to different models, e.g., in the context of quantum chemistry [3,4], many-body physics [5,6], or conformal bootstrap [34], often before or unaware of the general mathematical characterisation of non-commutative polynomial optimisation presented in [7,8].As it happens for variational methods, relaxations also provide values for other observables beyond energy, o R , but again with no control about whether they are close to or bound in any way the value in the ground state.Therefore, when combined, all what relaxations and variational methods define is an energy interval in which the searched groundstate energy lies.
Before concluding this part, it is worth mentioning that SDP relaxations of polynomial optimisation problems have a long tradition in the context of classical spin systems.In this case, the mathematical framework is the one of commutative polynomial optimisation, such that: where p and g j are again polynomials but now over (classical, namely commuting) variables x = (x 1 , . . ., x n ) and one deals with standard positivity constraints.Many classical spin systems have the form of polynomials over spin variables σ i such that σ 2 i = 1, again a polynomial constraint.A hierarchy of SDP relaxations of commutative polynomial optimisation problems was introduced in [35] and, in fact, the formalism of [8] used in this work can be understood as the extension of the construction in [35] to the non-commutative case, see [8].

III. CERTIFICATION OF GROUND-STATE PROPERTIES
Variational methods and relaxations of polynomial problems are often seen as two complementary approaches that allow one to bound the ground-state energy from above and below.The main point of our work is to show that their combination is much richer than expected, as together they can be used to derive certifiable bounds on any observable of interest in the ground state.
The idea is quite simple and was also discussed in [11].As mentioned above, the ground-state energy problem can be seen as an instance of polynomial optimisation because the Hamiltonian can be expressed as polynomials of some operators X i .For instance, for finite dimensional systems, it is enough to take as X i a basis for the space of matrices at each site, say Pauli matrices for qubit systems.In fact, any observable of interest O can be expressed as polynomial on these operators and bounds on it can be derived through the NPA formalism.Now, to restrict this optimisation to a region close to the ground state, one can use the best upper bound E A to the ground-state energy derived through variational methods, as well as the best lower bound E R derived through relaxations.The resulting optimisation reads: The different SDP relaxations of this minimisation provide a sequence of lower bounds to the actual value of the observable in the ground state, o (1) Note that in this case, the asymptotic value o (∞) is only guaranteed to be equal to the actual ground-state value o GS if E A = E GS .Finally, upper bounds o UB can be derived just by replacing the minimisation in ( 6) by a maximisation, obtaining the announced certifiable bounds for any observable in the ground state, To illustrate the power of this method, we apply it in what follows to several paradigmatic Heisenberg models for spin-1/2 systems.

IV. APPLICATIONS AND RESULTS
We present several implementations of the method to obtain certified bounds on ground-state observables for various Heisenberg models in one and two spatial dimensions [36].Generic Heisenberg models are defined by Hamiltonians of the form: where i ∈ {1, 2, . . ., N } label the lattice sites, while the couplings J ij implicitly define the lattice geometry and σ a i are the Pauli matrices acting on site i.The 1/4 prefactor follows standard condensed-matter conventions, where Hamiltonians are typically defined in terms of spin operators s a i = σ a i /2 instead of Pauli matrices.We shall consider four different geometries: 1.The Heisenberg model with first-neighbour interactions on a 1D lattice, J ij = δ j,i+1 , with periodic boundary conditions (PBC), namely we use the convention that N + 1 ≡ 1 for the i, j labels.
2. A 1D lattice with first-and second-neighbour couplings, , where the J 2 term induces geometric frustration.
Here, lattice sites are labelled by i = (x, y) with x, y ∈ {1, 2, . . ., L} (so that N = L 2 ), and couplings are of the form . We take PBC, namely L + 1 ≡ 1 for both x and y labels.
4. A 2D square lattice with first-and secondneighbour (frustration-inducing) couplings, where second-neighbours are along the diagonal of elementary square plaquettes, namely, extra couplings of the form In all cases, we obtain certified lower bounds on the ground state energy, as well as upper and lower bounds on relevant observables in the ground state, typically on spin-spin correlation functions.

A. Algorithmic considerations
We begin with discussing the concrete SDP algorihm tailored to generic Heisenberg models [Eq.(7)].In particular, we briefly discuss how to reduce the size of SDP relaxations by exploiting algebraic structures of the model, which is crucial in order to obtain the optimal results given some computational resource.More details are given in Appendix A.
The ground state energy of the Heisenberg model is the optimum of the following non-commutative polynomial optimisation problem: SDP relaxations can then be applied to Eq. ( 8) to obtain lower bounds of the ground-state energy, by replacing the optimisation over quantum many-body states |ψ⟩ by an optimisation over moment matrices.Specifically, suppose that B = {u m } is a monomial list, i.e., a subset of monomials of the form u m = σ a1 i1 . . .σ ap m ip m with respect to the non-commuting variables {σ a i } i=1,...,N ;a∈{x,y,z} , where p m is the degree of the monomial.Suppose also that the energy ⟨H⟩ can be expressed as a linear combination of the entries of the moment matrix M indexed by B with [M] uv = ⟨u † v⟩.Then a SDP relaxation to Eq. ( 8) is given by: min M obeys some moment replacement rules.(9) In particular, the equality constraints in Eq. ( 8) give rise to corresponding replacement rules on monomials, allowing one to reduce them to the normal form NF(u It follows that the moment matrix M satisfies the moment replacement rule: ⟨u⟩ = ⟨NF(u)⟩ for all entries u of M. Note that ( 9) is a complex SDP.To reformulate (9) as a SDP over real numbers, we refer the reader to [37].
Further symmetry considerations on the concrete considered models allow one to drastically reduce the number of independent variables in the moment matrix.For instance, given the symmetry of Heisenberg models under global rotations of the spins, correlations in the ground state are of the form (1/4)⟨σ a i σ b j ⟩ = δ a,b C ij .Another relevant symmetry is translation invariance which implies that correlation functions only depend on the relative position of the spins.Details on the technical implementation of those and other symmetries to reduce the computational complexity of the SDP algorithm are given in Appendix A.

B. Heisenberg chain
The Heisenberg chain is defined by the Hamiltonian: where N + 1 ≡ 1 in order to implement PBC.The ground state is critical (namely: gapless in the thermodynamic limit) and displays antiferromagnetic correlations decaying as a power-law with distance, ⟨s a i s a i+r ⟩ = (1/4)⟨σ a i σ a i+r ⟩ = C r ∼ (−1) r /r α with some exponent α [38].
Ground-state energy.-Theground-state energy per spin is given by e PBC (N ) = ⟨H⟩/N = 3C 1 .In Fig. 1, we plot the best lower bound of e PBC as obtained by our SDP relaxation for up to N = 100 spins.As a comparison, we plot the (quasi-)exact energy as obtained by DMRG simulations.It is also of interest to compare our SDP lower bound with the Anderson bound, which is obtained by exact diagonalisation on a system with open boundary conditions (OBC), e OBC (L) ≤ e PBC (N ) for all N > L (see Appendix B for details on the Anderson bound).In Fig. 1, the SDP bound is seen to vastly outperform the Anderson bound (which is in fact estimated with DMRG for the sake of illustration for up to N = 100, as beyond a few tens of qubits exact diagonaisation is out of reach).
To build the moment matrix in the SDP construction, we use all monomials of the form: with j ∈ {1, 2, . . ., r} and a, b, c, d ∈ {x, y, z} (all different monomials appearing only once).For each size N , we have chosen r as large as possible compatible with memory limitations, namely r = N 2 for N ≤ 60, and r = 20 for N = 80, 100.Furthermore for N = 100 we discard all degree-four monomials σ a i σ b i+1 σ c i+2 σ d i+3 in order to allow for more degree-two monomials.For the sake of completeness, the data plotted in Fig. 1 are also reported in Table III.Combining both the DMRG upper bound e DMRG and the SDP lower bound e SDP allows us to sandwich the exact ground-state energy with a relative accuracy that remains below 10 −3 up to N = 100 spins.In contrast, previous works have achieved no better than a few percents accuracy for comparable system sizes [5,6,39].The small energy gap between the DMRG (variational upper bound) and the SDP (certified lower bound) therefore certifies both the expected good performance of DMRG to approximate the actual 1D ground state and, in turn, also the good performance of the implemented SDP relaxation.III).Upper and lower bounds are derived, respectively, through DMRG and the implemented polynomial relaxation.For comparison, we show the expected Anderson bound, that is, the lower bound one would obtain by exactly solving the same system with OBC.This estimation is also computed through DMRG, so that it is possible to plot system sizes that are out of reach for exact diagonalisation.

C. Heisenberg chain with second-neighbour couplings
Our second application is the Heisenberg chain including both first-and second-neighbour couplings, namely the so-called J 1 − J 2 Heisenberg model: with PBC (in our convention the first-neighbour coupling is J 1 = 1).The J 2 term induces geometric frustration, leading to the sign problem in quantum Monte Carlo methods and to a richer phase diagram.The model was investigated in early days of DMRG simulations [40], and represents a cornerstone in the study of quantum magnetism, motivating the development of various variational wave-functions.In particular, it is predicted that for J 2 < J 2,c = 0.241167 . . ., the spin correlation length is infinite, and correlations decay as a power-law as in the J 2 = 0 limit [40].For J 2 > J 2,c , a gap opens and the system spontaneously forms dimers among first neighbours.In particular, at J 2 = 0.5, the two exact ground states are products of Bell pairs among first neighbours [41].For larger values of J 2 , more complex correlation patterns emerge, with both long-range dimer-dimer correlations and finite-range spiral spin correlations [40].Those predictions are based on DMRG (hence, variational) numerical simulations.Here, in contrast, we investigate the ability of SDP techniques to offer relevant lower bounds to the ground-state energy, as well as certified bounds on spin correlations in the ground state -something that, to our knowledge, no other approach can provide.In particular, we certify a change of sign for the second-neighbour spin correlations for J 2 > 0.5 (see Fig. 4).
Ground state energy.-InFig. 2, we plot the best lower bound of the ground-state energy for a system of size N = 40, as compared to the DMRG value (the data are reported in Table IV; see also Table V for the lower bounds computed for N = 100).Our compromise for the choice of monomials is different for small and large values of J 2 .For J 2 ≤ 1, the monomials are the same as for J 2 = 0, namely: with j ∈ {1, 2, . . ., r} and a, b, c, d ∈ {x, y, z}.For J 2 > 1, to better capture the effect of frustration, our (heuristic yet efficient) choice is: with j ∈ {1, 2, . . ., r} and a, b, c, d ∈ {x, y, z}.
As can be seen in Table IV, for J 2 ≲ 0.5 we obtain a relative accuracy of 10 −3 , and for all values of J 2 the relative accuracy remains better than 0.016.As  IV).Upper and lower bounds are derived, respectively, through DMRG and the implemented polynomial relaxation.Inset: the relative accuracy remains better than 0.016 for all values of J 2 .
expected, the largest gap appears at J 2 = 1.0,where the two couplings become comparable and there is competition between them.
Individual terms in the Hamiltonian.-Asmentioned, the SDP approach allows one to obtain certified bounds on relevant observables in the ground-state beyond the energy.In order to do so, we constrain the energy to lie in-between the DMRG upper bound and the SDP lowerbound.As a first application, we compute bounds on the first-neighbour spin correlations C 1 (see Fig. 3), as well as the second-neighbour spin correlation C 2 (see Fig. 4), namely both individual terms composing the Hamiltonian.For the sake of comparison, we also plot the results obtained through DMRG calculations, which are expected to be very close to the exact value.The derived lower and upper bounds certify that: • First-neighbour correlations remain antiferromagnetic, C 1 < 0, for all values of J 2 , as its upper bound is always negative (see Fig. 3 and Table VI); • at J 2 = 0.5, the second-neighbour correlations change from ferromagnetic (C 2 > 0) to antiferromagnetic (C 2 < 0) (see Fig. 4 and Table VII).This is a non-trivial qualitative information, illustrating the competition between the J 1 term which favors staggered correlations among first neighbours (namely C 1 < 0 and C 2 > 0) and the J 2 term which favors staggered correlations among second neighbours (namely C 2 < 0).
These findings are fully compatible with the DMRG results, but recall that the latter cannot provide any certification about these properties.This is our first illustration of how physically relevant correlation properties in the ground state can be certified using SDP relaxations.VII).

SDP Lower
Spin correlations.-Wethen study the ability of the SDP approach to bound the spin correlation function at larger distance.In particular, as mentioned, for values of J 2 < J 2,c , one expects that the system develops antiferromagnetic (that is, staggered) spin correlations which decay as a power-law with distance, as in the J 2 = 0 limit [40].In order to explore the potentiality of SDP relaxations to capture such quasi-long-range order the ground state, we compute bounds on the spin-spin correlations as a function of distance for a fixed system size of N = 40.We consider both J 2 = 0.2 < J 2,c (Fig. 5 and Table VIII) and J 2 = 1.0 > J 2,c (Fig. 6 and Table IX).
As can be seen, the SDP upper and lower bounds tightly sandwich the DMRG value at small distances, while they  IX).Note that the bounds at i = 19 show a remarkable and surprising improvement.We have however verified that the obtained SDP solutions define feasible points and hence provide valid lower and upper bounds.
become looser at larger distances.Yet, one sees that • For J 2 = 0.2, the SDP bounds are tight enough to certify the staggered sign structure of the correlation function up to the maximal distance i = N/2 (see Fig. 5 and Table VIII).Indeed, both the lower and upper bounds change sign with the distance.
• For J 1 = 1.0,SDP bounds certify instead a qualitatively different spatial structure of spin correlations at short distance, while they become much looser at larger distance.
It is important to remark that we have not attempted here to optimise the choice of monomials to best capture the correlation function at large distance; instead we have kept the same monomials as for tightly bounding the energy, which especially constrain correlations at short distances.One can expect to get tighter bounds by tailoring the monomial list to the observable to be certified.We come back to this point below.Carlo (data in Table X).
We now move to the most challenging case of twodimensional systems.As above, we start with the Heisenberg model, but now on a square lattice: where (i, j) label the position of the spins on a square lattice with PBC (L + 1 ≡ 1).In contrast to the 1D model, it is expected that the square-lattice Heisenberg model spontaneously breaks the SU (2) symmetry in the thermodynamic limit and displays true long-range antiferromagnetic order in the ground state, implying in particular that C(L/2, L/2) → cst.> 0 for L → ∞.As further discussed below, while we do certify this property for L ≤ 8, we cannot reliably extrapolate the obtained SDP bounds to the thermodynamic limit.
We compare the derived bounds to the quantum Monte Carlo data from [42], which are expected to be equal to the exact values (up to statistical error bars which are negligible on the scale of our comparison) (see Fig. 7 and Table X).The energy gaps are now larger than in the one-dimensionl case, but the relative accuracy of the SDP lower bound of the energy is still about 0.01 as compared to the quantum Monte Carlo result.Long-range order.-Wenow focus on other groundstate properties beyond energy and, in particular, on long-range correlations.In order to investigate the possibility to certify spontaneous symmetry breaking and the associated long-range antiferromagnetic order in the ground state, we compute bounds on the correlation at maximal distance C(L/2, L/2) (see Fig. 8 and Table XI).For all the computed sizes, the lower bound on C(L/2, L/2) remains positive, hence certifying the presence of long-range order on those sizes.However, similarly to the case of the J 1 − J 2 model at J 2 = 0.2 (Section IV C), for increasing system size the SDP bounds on C(L/2, L/2) become increasingly looser.It is therefore not possible to argue that C(L/2, L/2) will remain positive for L > 8 (corresponding to 1/L < 0.125 on the figure) from scaling arguments on the derived lower bounds.
Note again that, to derive these bounds on long-range correlations, we have used the same monomial list as for optimising the ground state energy.In future works one may instead tailor the choice of the monomials to better bound the correlation function at large distance, which can be expected to offer some improvement.
E. J1 − J2 square lattice Heisenberg model As a last example, we consider the J 1 − J 2 Heisenberg model on a square lattice: +J 2 (σ a (i+1,j+1) + σ a (i+1,j−1) ) , (13) with PBC.The J 2 terms favors antiferromagnetic correlations along the diagonals of the square lattice, which are incompatible with the correlations favored by the first-neighbour J 1 = 1 term and leads to frustration.
As for the J 1 − J 2 Heisenberg chain, this model is not amenable to quantum Monte Carlo due to the sign problem.Several variational methods based on Ansatz wavefunctions have however been applied to this paradigmatic model of frustrated quantum magnetism, sometimes obtaining conflicting results due to a complex energy landscape with various ground-state candidates which are close in energy yet with incompatible forms of order [43][44][45][46][47][48].
Ground-state energy.-Again,we first compute SDP lower bounds on the energy, which complement variational methods.We present results for L = 6, 8 in the Appendix G 1 (see Figures 12 and 13), and for L = 10 in Figure 9 (data are respectively given in Tables XII, XIII and XIV).Notice in particular that size L = 10 (namely, N = 100 qubits) is not achievable with exact methods, so that combining upper-and lower bounds become very relevant to constrain ground-state properties.Combining variational upper bounds and SDP lower bounds allows us to sandwich the true ground-state energy with a few percent of relative accuracy (Table XIV Spin correlations.-Finally,the SDP approach can also be applied to deliver certified bounds on relevant observables in regimes inaccessible to exact numerical methods, such as the J 1 − J 2 model on a square lattice for L = 10.Constraining the energy to lie in-between the maximal lower bound (as obtained by the SDP) and the minimal upper bound (as obtained using variational methods), we obtain certified bounds on first-and second-neighbour (diagonal) correlations in the exact ground state.The monomial list that we use to bound both the energy and spin correlations is the same as the one in Section IV D for the square lattice Heisenberg model.
The results are respectively displayed in Fig. 10 and Fig. 11 (data in Tables XV and XVI).We emphasize that this frustrated model for N = 100 spins is well beyond the capabilities of known exact methods such as exact diagonalisation, so that the certified bounds offered by the SDP approach are especially insightful.We notice in particular that SDP bounds are sufficiently accurate to certify: • Correlations C(0, 1) remain antiferromagnetic for all studied values of J 2 , as the computed upper bound is always negative.
• Second-neighbour correlations experience a change of sign while varying the J 2 /J 1 ratio, a behavior reminiscent of the 1D model studied in Section IV C. The transition occurs for a value of J 2 in the range (0.45, 0.6).
Again, this type of certification is impossible with previous approaches.XV).

V. DISCUSSIONS
In this work we have shown how SDP relaxations of polynomial optimisation problems when combined with upper bounds obtained through variational methods can provide certifiable bounds on ground-state properties beyond energy.We have illustrated the potentialities of the method in 1D and 2D Heisenberg models.The choice is motivated by their rich phenomenology, the existence of previous results to benchmark our results, and their symmetries, which allow us reaching large system sizes.However, the method is general and can be applied to essentially any many-body Hamiltonian.In fact, a first natural continuation of our results is to apply the introduced approach to other relevant models in physics, for instance for fermions.Note also that while our work has focused on finite-size systems, it is also possible to adapt it to systems in the thermodynamic limit, following similar ideas as discussed in [49].
Another direction that deserves further considerations is to understand how to use the knowledge of the model under study to improve the obtained bounds.Our work already does this by exploiting the symmetries and sparsity of the considered models to significantly improve the scalability.But one can also study how to combine our approach with different relaxation techniques, such as those of [49].Another way of improving our results is by adapting the choice of monomials to the observable of interest.For instance, in our work, we choose the monomials based on the local operators appearing in the Hamiltonian.However, once the bounds on the energy are obtained, one could have modified the monomials in the SDP when studying long-range correlations.We leave the study of this possibility for further work.Note also that machine learning can also be employed to choose the monomials, as done in [50].
On the numerical side, some of the lower bounds based on a solution returned by a SDP solver may come together with unsatisfying numerical feasibility status.Therefore, another interesting research direction is to obtain truly certified lower bounds based on exact rational arithmetic.For this, one could design a post-processing method, relying either on rigorous interval arithmetic or rounding-projection techniques, in the same spirit as in [51,52].Also, the SDPs arising from the NPA hierarchy possess some special structures (e.g., low-rank opti-mal solutions, unit diagonal) which could be exploited to design more efficient SDP algorithms as in [53].We can thus rely on a structure-exploiting SDP solver to approach models of larger size in the future.
To conclude, our work demonstrates how to derive bounds on ground-state properties of quantum manybody Hamiltonians beyond energy.In doing so, it opens many different perspectives for future research.Our vision is that the considered techniques will become a standard tool to complement and certify result of variational methods, so far the dominant tool to study quantum many-body problems.

VI. ACKNOWLEDGEMENTS
This work is supported by the ERC AdG CERQUTE, the Government of Spain (NextGenerationEU PRTR-C17.I1 and Quantum in Spain, Severo Ochoa CEX2019-000910-S), Fundació Cellex, Fundació Mir-Puig, Generalitat de Catalunya (CERCA programme), the AXA Chair in Quantum Information Science, EU Quantera project Veriqtas, the National Natural Science Foundation of China under Grant No. 12201618, the European Union Horizon's 2020 research and innovation programme under the Marie Sklodowska Curie grant agreement No 101031549 (QuoMoDys).Research at the Perimeter Institute for Theoretical Physics is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Research, Innovation and Science.VM was supported by the FastQI grant funded by the Institut Quantique Occitan, the PHC Proteus grant 46195TA, the European Union's Horizon 2020 research and innovation programme under the Marie Sk lodowska-Curie Actions, grant agreement 813211 (POEMA), by the AI Interdisciplinary Institute ANITI funding, through the French "Investing for the Future PIA3" program under the Grant agreement n • ANR-19-PI3A-0004 as well as by the National Research Foundation, Prime Minister's Office, Singapore under its Campus for Research Excellence and Technological Enterprise (CREATE) programme.
Appendix A: SDP reductions by exploiting structure In this appendix, we further discuss how to implement symmetries of the problem in the SDP algorithm to reduce the number of free variables in the implementation.We start be restating the main problem for completeness.
The ground state energy of the Heisenberg model is the optimum of the following non-commutative polynomial optimisation problem: (A1) The NPA hierarchy [8] can be then applied to (A1) yielding a non-decreasing sequence of lower bounds on the ground state energy.Specifically, suppose that B d is a monomial basis (i.e., a subset of monomials w.r.t. the non-commutating variables {σ a i } i=1,...,N,a∈{x,y,z} ) up to degree d.Then the d-th order moment relaxation of the NPA hierarchy for (A1) is given by min M d obeys some moment replacement rules, (A2) where M d is the moment matrix indexed by B d with [M d ] uv = ⟨u † v⟩.Note that the equality constraints in (A1) give rise to the following replacement rules on monomials: For any monomial u, by applying the above replacement rules, we can reduce it to the normal form NF(u It follows that the moment matrix M d satisfies the moment replacement rule: ⟨u⟩ = ⟨NF(u)⟩ for all entries u of M d .

Sparsity
In order to exploit the sparsity of the Heisenberg model, for each degree d, we pick monomials that are supported on contiguous sites: Then at relaxation order d, we use the sparse monomial basis B d = ∪ d i=0 P i instead of the full monomial basis.Moreover, to capture long-range correlations, we also include the monomials of form σ a i σ b i+j with j = 2, . . ., r and a, b ∈ {x, y, z} in the monomial basis.The resulting SDP relaxation is more efficient to solve but possibly leads to more conservative lower bounds.We emphasize that similar reduction of the monomial basis already appeared in the related literature.In the context of quantum information theory, we refer to [54] where the authors obtain upper bounds on maximal violations of Bell inequalities after random selection of a subset of monomials with given degrees.For general (non-)commutative polynomial optimisation problems, one can exploit either correlative sparsity [55,56], occurring when there are few correlations between the variables of the input problem, or term sparsity [57,58], occurring when there are a small number of terms involved in the input problem by comparison with the fully dense case.The interested reader is referred to [59] for a recent monograph on this topic.

Symmetry a. Sign symmetry of the model
We can observe that the feasible set of (A1) is invariant under the substitution of two of the three variables, e.g., σ x i and σ y i of a given site into their opposite, e.g., −σ x i and −σ y i .In order for any objective functions of the form (7) to also be invariant, we need to consider the same substitutions for all the sites.There are therefore three substitutions: Note that s zx is the composition of s xy and s yz so we only need to consider the invariance under two of the three substitutions.
For each monomial m, s xy (m) (resp.s yz (m)) is either m or −m.Similarly to [60, Section III.C], for each monomial m, we consider its signature as the vector (s xy (m)/m, s yz (m)/m) ∈ {−1, 1} 2 .Consider that the list of monomials is such that each of the four groups of monomials with the same signature appears contiguously.Notice that a product of two monomials is invariant under the sign symmetries if and only if is the product of two monomials of the same signature.The moments of symmetric monomials therefore form a block diagonal structure of 4 blocks in the moment matrix.We can reduce the large positive semidefinite matrix into 4 smaller positive semidefinite matrix as shown in [60,Theorem 4] in the commutative case.
For example, when d = 1, the partition is given by Table I.When d = 2, the partition is given by Table II.The Hamiltonian H in (A1) may have more sign symmetries.For example, consider the Heisenberg model, the Hamiltonian H is invariant under the substitutions: Besides the zero entries given in Section A 2 a, these additional sign symmetries of the Hamiltonian yield extra zero entries of the moment matrix: ⟨u⟩ = 0 if NF(u) is variant under the transformations (A5).
Note that symmetry reduction usually applies for any convex problem with symmetric objective and symmetric feasible set.In this case, the symmetry of the feasible set is not apparent because the replacement rules are not symmetric, e.g., replacing σ x i σ y i by iσ z i is not symmetric under any of the substitutions (A5).We can show the symmetry of the set of sums of Hermitian squares using (1) their connection with the set of strictly positive Hermitian elements over the quotient ring (detailed below) and ( 2) the fact that the Hamiltonian is invariant under the substitutions (A5).
Indeed, let us denote the non-commutative polynomial ring by C⟨{σ x i , σ y i , σ z i } N i=1 ⟩ and the ideal generated by the equality constraints of (A1) by I. Let Ω := {NF(u) | u is a monomial in {σ x i , σ y i , σ z i } N i=1 }.Then the optimization problem (A1) is equivalent to the unconstrained optimization problem: min {|ψ⟩,σ a i } ⟨ψ|H|ψ⟩, considered in the quotient ring C⟨{σ x i , σ y i , σ z i } N i=1 ⟩/I ∼ = C⟨Ω⟩.Let us denote by S the group of additional sign symmetries given by (A5) and by Σ S the set of sums of Hermitian squares of C⟨Ω⟩, that are invariant under S after conversion to normal form, i.e., elements of the form h = j p † j p j , p j ∈ C⟨Ω⟩, such that s(NF(h)) = NF(h) for any s ∈ S. Then one can show that any strictly positive Hermitian element of C⟨Ω⟩ that is invariant under S lies in Σ S by [9], the proof being very similar to the one of [61, Proposition 3.1].Here "strict positivity" should be understood as strict positivity over all possible evaluations in Hilbert spaces, as detailed, e.g., in [61,Section 2.2].By duality one considers optimization over linear functionals nonnegative on Σ S , which leads to an equivalent formulation of the above unconstrained optimization problem: From any ℓ feasible for the above problem (A6), let us define the linear functional ℓ S : C⟨Ω⟩ → C by ℓ S (p) = 1/|S| s∈S ℓ(s(p)) for all p ∈ C⟨Ω⟩.Then it is clear that ℓ S (H) = ℓ(H) since H is invariant under the action of S. In addition, ℓ S (1) = 1, ℓ S (p † ) = ℓ(p) * for all p ∈ C⟨Ω⟩ and ℓ S (q) ≥ 0 for all q ∈ Σ S .To prove the latter fact, we used the fact that one has s(q) = q for any s ∈ S and any q ∈ Σ S written in normal form.Overall this shows that ℓ S is feasible for (A6) and yields the same objective value than the one with ℓ.
To conclude, at each relaxation (A2) one can restrict ourselves to optimizing over linear functionals vanishing on variant elements of C⟨Ω⟩ under the transformations (A5).This boils down to setting every entry of the Hermitian moment matrix M = [ℓ(v † w)] v,w∈Ω from (A2) to ⟨u⟩ = 0 if NF(u) is variant under the transformations (A5).

c. Translation symmetry
The translation symmetry of (A1) comes from that the Hamiltonian H is invariant under any translation of sites, which implies where υ : i −→ i+r denotes a translation of sites with r ∈ {1, . . ., N }.This together with the PBC imposes a block structure on the moment matrix M d where each block is an Hermitian circulant matrix as long as the monomial  as SDP data for N = 100 spins (Table V).In both tables, the last column indicates the degree d of the monomials to construct the moment matrix, as explained in Appendix A. We also provide numerical data for the spin-spin correlations in a system with N = 40 spins with PBC.
We then provide the numerical data for the spin-spin correlation as a function of distance, for both J 2 = 0.2       In Table XI we provide data on the spin correlation C(L/2, L/2), namely at maximal distance along the diagonal of the square lattice (main Fig. 8).It is expected that this correlation remains nonzero in the thermodynamic limit, corresponding to antiferromagnetic longrange order in the ground state.In this appendix, we provide numeral data regarding the J 1 − J 2 Heisenberg model on a square lattice (size N = L × L and PBC), as discussed in Section IV E in the main text.

Ground-state energy
We first present data for the ground-state energy, as compared with various variational methods employed in previous works in the literature.
In Table XII we present data for L = 6.We compare the SDP lower-bound with neural-network (NN) Ansatz wavefunctions [48] and exact results [63].The same data are plotted in Fig. 12.
In Table XIII     Heisenberg model (L = 10).Last column: relative difference between the best variational upper bound and the SDP lower bound.

First-and second-neighbour spin-spin correlations
We finally provide SDP bounds for first-and secondneighbour spin-spin correlations as a function of the second-neighbour coupling J 2 , for a system of size L = 10.Data are respectively displayed in Table XV physical observables O N in the ground state: |ψ (N ) GS ⟩, see also Eq. (1), namely o (N ) GS = ⟨ψ (N ) GS |O N |ψ (N ) GS ⟩.Often, one is also interested in the thermodynamic limit of an infinite number of particles, which are typically inferred using scaling considerations by studying the dependence of o (N ) GS with N .In what follows, we often remove the dependence on N to simplify the notation.

FIG. 1 :
FIG.1: Ground-state energy per particle in the Heisenberg chain (data in TableIII).Upper and lower bounds are derived, respectively, through DMRG and the implemented polynomial relaxation.For comparison, we show the expected Anderson bound, that is, the lower bound one would obtain by exactly solving the same system with OBC.This estimation is also computed through DMRG, so that it is possible to plot system sizes that are out of reach for exact diagonalisation.

FIG. 2 :
FIG.2: Ground-state energy per particle in the J 1 − J 2Heisenberg chain (N = 40; data in TableIV).Upper and lower bounds are derived, respectively, through DMRG and the implemented polynomial relaxation.Inset: the relative accuracy remains better than 0.016 for all values of J 2 .

FIG. 5 :FIG. 6 :
FIG.5: Spin-spin correlation in the J 1 − J 2 Heisenberg chain for J 2 = 0.2 and system size N = 40 (data in tableVIII).The staggered sign structure is certified at all distances.

FIG. 7 :
FIG.7: Ground-state energy in the square lattice Heisenberg model, and comparison to quantum Monte Carlo (data in TableX).
(corresponding to Fig. 10 in the main text) and Table XVI (corresponding to Fig. 11 in the main text).
that is, the set of expectation values {⟨h i ⟩} that can be obtained from a quantum state |ψ⟩.In particular, when the system is translationally-invariant, all expectations values ⟨h i ⟩ are identical, and the minimisation is over a single real number: the expectation value of the operator h i .As such, the problem may look simple, but the characterization of M Q is generically very hard.For instance, for the case of Hamiltonians with local interactions, that is, when the operators h i act nontrivially on a reduced set of particles indexed by i, the expectation values read ⟨h i ⟩ = tr(h i ρ i ), where {ρ i } i is a set of reduced states acting on the different subsystems indexed by i that are compatible with a global N -particle state |ψ⟩. ).

TABLE I :
Table of monomials indexing the sign symmetry blocks when d = 1.

TABLE II :
Table of monomials indexing the sign symmetry blocks when d = 2.
b. Sign symmetry of the Hamiltonian

TABLE III :
Heisenberg chain energy in a system of N spins.r denotes the maximal distance between spins for two-body terms in the monomial list (see main text, section IV B).Note that the relative accuracy (third column) is below 10 −3 for all sizes.

TABLE IV :
Ground-state energy in the Heisenberg chain with second-neighbour (J 2 ) for N = 40 spins, as evaluated by SDP and DMRG methods.

TABLE V :
Ground-state energy lower-bound in the Heisenberg chain with second-neighbour (J 2 ) for N = 100 spins, as evaluated by the SDP algorithm.

TABLE VIII :
Lower and upper SDP bounds for the spin-spin correlator at distance i in the Heisenberg chain with second-neighbour couplings (J 2 = 0.2 and size N = 40).

TABLE IX :
Lower and upper SDP bounds for the spin-spin correlator at distance i in the Heisenberg chain with second-neighbour couplings (J 2 = 1.0 and size N = 40).

TABLE X :
SDP lower bound on the energy of the square lattice Heisenberg model as compared to quantum Monte Carlo results.

TABLE XI :
SDP lower and upper bound for the spin correlations at maximum distance in the square lattice Heisenberg model, sandwiching quantum Monte Carlo results.
Appendix G: Square-lattice J1 − J2 Heisenberg model we present data for L = 8.We compare the SDP lower-bound with variational Monte Carlo on

TABLE XII :
[47]nd-state energy for the square-lattice J 1 − J 2 Heisenberg model (L = 6).Last column: relative difference between the best variational upper bound and the SDP lower bound.FIG.12:Energylowerbounds for the 2D J 1 − J 2 Heisenberg model on a square with L = 6 (data in TableXII).Ansatz wavefunctions (VMC) and DMRG computations[47].The same data are plotted in Fig.13.

TABLE XIII :
Energy lower bounds for the square-lattice J 1 − J 2 Heisenberg model (L = 8).Last column: relative difference between the best variational upper bound and the SDP lower bound.In Table XIV we provide numerical data for Fig.9in the main text, namely ground-state energy for L = 10 for both SDP, NN and DMRG methods.FIG.13:Energylowerbounds for the square-lattice J 1 − J 2 Heisenberg model with L = 8 (data in TableXIII).

TABLE XIV :
Energy lower bound (SDP) and upper bounds (NN and DMRG) for the square-lattice J 1 − J 2