Energy as a detector of nonlocality of many-body spin systems

We present a method to show that low-energy states of quantum many-body interacting systems in one spatial dimension are nonlocal. We assign a Bell inequality to the Hamiltonian of the system in a natural way and we efficiently find its classical bound using dynamic programming. The Bell inequality is such that its quantum value for a given state, and for appropriate observables, corresponds to the energy of the state. Thus, the presence of nonlocal correlations can be certified for states of low enough energy. The method can also be used to optimize certain Bell inequalities: in the translationally invariant (TI) case, we provide an exponentially faster computation of the classical bound and analytically closed expressions of the quantum value for appropriate observables and Hamiltonians. The power and generality of our method is illustrated through four representative examples: a tight TI inequality for 8 parties, a quasi TI uniparametric inequality for any even number of parties, ground states of spin-glass systems, and a non-integrable interacting XXZ-like Hamiltonian. Our work opens the possibility for the use of low-energy states of commonly studied Hamiltonians as multipartite resources for quantum information protocols that require nonlocality.


I. INTRODUCTION
Nonlocality is a fundamental property of Nature in which the statistics obtained by performing some local measurements on some composite quantum systems cannot be reproduced by any local hidden variable model [1]. These so-called nonlocal correlations cannot be mimicked by any local deterministic strategy, even if assisted by shared randomness [2]. Nonlocality is detected by the violation of a Bell inequality [3] and it has recently been demonstrated in three loopholefree Bell experiments [4][5][6]. Detection of nonlocality is a sufficient condition to demonstrate, in a device-independent (DI) way, that the state producing such correlations is entangled. From an operational point of view, nonlocality is a resource that enables the implementation of DI quantum information protocols, such as DI quantum key distribution [7,8], DI randomness expansion [9] or amplification [10,11], or DI selftesting [12,13].
The study of quantum many-body systems has benefited during the last decades from insights of the field of quantum information; in particular concerning the understanding of their correlations [14,15]. This has, however, mostly focused on the study and experimental detection of entanglement [16][17][18], while the role of nonlocal correlations, which are stronger, remains rather unexplored. There are at least three reasons for that: First, the known Bell inequalities for multipartite systems involve correlations among many particles [19][20][21][22][23], thus rendering their measurement a formidably challenging task. Second, the mathematical characterization of nonlocal correlations is an NP-hard problem [24]. Third, the size of the description of multipartite quantum states grows, in general, exponentially with the system size, posing a strong barrier to the analysis of the quantum correlations in large systems. However, recent advances [25,26] have shown that, by measuring only one-and two-body correlation functions, nonlocality can be revealed in some multipartite quantum systems, opening the way to its detection in many-body systems [27] (see also [28]).
In this work we show that the ground states of some quantum spin Hamiltonians in one spatial dimension are nonlocal. We assign a Bell inequality to the given Hamiltonian in a natural way and we calculate its classical bound using dynamic programming. The Bell inequality is constructed in such a way that, for appropriate quantum observables, the Bell operator coincides with the Hamiltonian. The idea is that if the ground state energy is beyond the classical bound, this signals the presence of nonlocal correlations in the ground state. The ground state energy is computed by exact diagonalization using the Jordan-Wigner (JW) transformation, which maps the system of spins to a quadratic system of fermions. The method just presented can also be seen from the opposite point of view, namely, as a way to optimize certain classes of Bell inequalities for many-body systems under some quantum observables. We also study the translationally invariant (TI) setting, in which we provide an exponentially faster algorithm to find the classical bound and we obtain analytical results for the quantum value 1 . Then we illustrate our framework by applying it to four representative examples: First, a tight TI Bell inequality for 8 parties. Second, a quasi TI Bell inequality for any even number of parties. Third, we show that the ground state of an XY spin glass has nonlocal correlations in some parameter region. Finally, a non-integrable interacting XXZlike Hamiltonian to which we assign a variation of Gisin's Elegant Bell inequality [29] and we find its quantum value numerically using matrix product states and the density matrix renormalization group [30]. This shows that our method is not limited to Hamiltonians that can be mapped to a system 1 In this work, we refer to the quantum value as the expectation value of the ground state of the Bell operator under appropriate quantum observables (such that it matches the Hamiltonian). This should not be confused with the quantum bound of a Bell inequality, which is the infimum over all quantum states and observables arXiv:1607.06090v2 [quant-ph] 1 May 2017 of free fermions via the JW transformation, but it can be applied to any spin Hamiltonian with short-range interactions in one spatial dimension. The paper is organized as follows: In Section II we present our method. In Section III we show how to optimize certain classes of many-body Bell inequalities. In Section IV we particularize our results to the TI case. In Section V we illustrate our methods with four examples. Finally, in Section VI we conclude and explore future lines of research.

II. THE METHOD
In this section we present a method to analyze when the ground state of some spin Hamiltonians is nonlocal; namely, when the quantum value is beyond the classical bound. We first present the setting (Section II A). We then construct a Bell inequality from the given Hamiltonian (Section II B) and we compute its classical bound using dynamic programming (Section II B 1). Then we find the quantum value of the inequality (Section II C). To this end, we first review the exact diagonalization of quadratic fermionic Hamiltonians (Section II C 1) and then we relate the spin Hamiltonian and the fermionic Hamiltonian via the Jordan-Wigner (JW) transformation (Section II C 2).

A. The setting
We consider quantum spin-1/2 Hamiltonians of n particles in one spatial dimension (henceforth, one-dimensional) with periodic boundary conditions, with short-range interactions, up to R neighbors: where t (i) and t (i,r) α,β are real parameters, are the so-called String operators, σ (i) x , σ (i) y and σ (i) z are the Pauli Matrices acting on the i-th site and the indices of the sites are taken modulo n. We denote x + 1 := y for short.
On one hand, as we explain it in detail in Section II C, this choice of Hamiltonians is convenient from a purely mathematical perspective, as these/they can be exactly diagonalized via the JW transformation. On the other hand, Hamiltonians of the form (1) are general enough to include many cases of interest. For instance, the case R = 1 corresponds to a onedimensional spin-1/2 Hamiltonian with nearest-neighbors interactions, under a transverse magnetic field: Nevertheless, our method can also be applied to Hamiltonians with local interactions that do not rely on the String operator structure, at the price of having to compute their ground state energy numerically, as we illustrate it in Example V D.

B. Construction of a Bell inequality
In this section we study the classes of Bell inequalities that are relevant for our work. We construct Bell inequalities such that, for some quantum observables, the corresponding Bell operator B satisfies where β C ∈ R is the so-called classical bound and H is defined as in Eq. (1). We use the following convention in writing down Bell inequalities. We want to obtain Bell inequalities of the form I + β C ≥ 0. The part of the Bell operator that corresponds to I will be the Hamiltonian H, and thus, states with low enough energy, lower than −β C , will give a violation of the Bell inequality. The classical bound β C is defined as where the minimum is taken over all Local Hidden Variable Models (LHVM) (cf. Eq. [3]). Observe that the quantum state that minimizes the expectation value of B is the ground state of H. This motivates the study of Bell inequalities in the following scenario: We have n parties with m dichotomic observables with outcomes ±1 at their disposal. We denote the choices of measurements by k = (k 0 , . . . , k n−1 ), with 0 ≤ k i < m, and the outcomes they produce by a = (a 0 , . . . , a n−1 ) with a i = ±1. We denote by P (a|k) the vector of conditional probabilities collected when they perform the Bell experiment. Due to the no-signalling principle, the marginals observed by any subset of parties do not depend on the choices of measurements performed by the rest; thus P ({a i } i∈S |{k i } i∈S ) is well defined on any subset S. In the case of dichotomic measurements, one normally works with the correlators M where abusing notation we are now denoting k = (k i , . . . , k i+r ) and a = (a i , . . . , a i+r ).
If R > 1, we will consider m + 1 measurements, due to the σ z inbetween (cf. Eq. (2)). The Bell inequalities that are naturally tailored to a Bell operator of the form of Eq. (4) can be written as I + β C ≥ 0, where and γ (i) , γ (i,r) k,l are real parameters that depend on the t (i) and t (i,r) α,β of Eq. (1). Despite the fact that I contains up to (R + 1)body terms, its coefficients γ (i) and γ (i,r) k,l show that it is essentially a 2-body Bell inequality, since the measurement choice of the parties in the middle of the string is fixed to be m, in the sense that the number of coefficients γ (i) and γ (i,r) k,l is the same as in a 2-body Bell inequality. Note that the number of measurements m performed on the x-y plane will not affect the form of Eq. (1). The only measurement that is not performed in this plane is in the z direction; therefore, we treat it as a special case and we say we have m + 1 measurements.

The classical optimization
We now describe how to efficiently compute the classical bound β C of the Bell inequalities introduced in Eq. (7). It is well known that for a generic Bell inequality for n parties, m measurements and d outcomes the classical bound cannot be found efficiently, as it requires solving a linear program with d mn constraints [3]. The particular form of the Bell inequalities we are considering, however, allows us to find an algorithm that, in the many-body regime (i.e., for fixed d, m and R) has O(n) complexity. For simplicity, we consider a Bell inequality I of a slightly more general form than those of Eq. (7) and with d = 2 outcomes, where k = (k i , . . . , k i+r ) and 0 ≤ k j < m. After we have presented our method, it will become clear that there is no loss of generality in considering dichotomic measurements, as the result can be straightforwardly generalized to an arbitrary d.
To find β C , we need to optimize I over all local hidden variable models. By Fine's Theorem [2], it suffices to optimize I over all deterministic local strategies, in which the correlators M (i,r) k factorize as where M (i) ki can be ±1. Thus, where the minimum is taken over all possible assignments of M (i) ki to ±1 for all i and k. Let us first solve the case with Open Boundary Conditions (OBC); i.e., the case where γ (i,r) k = 0 when i + r ≥ n. We shall follow a dynamic programming procedure [31] that splits the minimization (10) into nested parts.
To this end, we represent a local deterministic strategy as a matrix M whose rows index the measurement choices and whose columns index the party, and the entry (k, i) is the value · · · · · · · · · . . . . . .
The representation of the matrix M at the step where we compute Ei for an inequality with dichotomic measurements with R = 2. The color (red or black) corresponds to the outcome assigned to that observable. At the i-th step, an optimal local deterministic strategy for parties 0 to i − 2 is already set, and we are minimizing the local deterministic strategy at the (i − 1)-th site (green squares). In order to perform the minimization, we need the values of the ith and (i + 1)-th parties to be fixed (blue triangles); i.e., we need to give a value to M (i,2) . We check the 2 m possible values for the observables at site i − 1 and we find an optimal assignment for them. assigned in the deterministic strategy to the k-th observable of the i-th party. Thus, M is a (m × n) matrix whose entries take integer values +1 or −1. Let M (i,R) denote the submatrix of M consisting of columns i to i + R − 1. The goal of the dynamic programming is to find an optimal M, which will give β C . This will be obtained recursively.
Let h i be the function defined for i > 0 as Note that because of Eq. (9), h i (M (i−1,R+1) ) is a real number. Now we define a recursive function E i which contains the optimization up to the (i − 1)-th site. Explicitly, E 0 (M (0,R) ) := 0 and for i > 0. Note that E i optimizes the local deterministic strategy on the (i − 1)-th party, which amounts to choosing the optimal values of the (i − 1)-th column of M. This naturally depends on the next R columns, which we need to consider as variables and calculate E i for all their possible values. Therefore, to efficiently evaluate E i (M (i,R) ), we only need to access the stored values of E i−1 on M (i−1,R) and not E i−2 and so on. The computation of E i (M (i,R) ) thus requires evaluation of the 2 m different deterministic local strategies corresponding to the values M (i−1) k (see Figure 1).
The classical bound (10) is obtained as the end product of this minimization procedure, namely Note that E n is actually independent of M (n,R) because we are in the OBC case, so all the γ k 's that would extend beyond the n-th party are zero. This procedure can be easily generalized if the Bell inequality has d > 2 outcomes. In this case, one has to take into account that the notion of correlators introduced in Eq. (6) is no longer well defined (as the one-to-one correspondence between probabilities P (a|k) and correlators M (i,r) k no longer holds). Thus, one needs to express the Bell inequality in terms of probabilities: where a = (a i , . . . , a i+r ) and k has the same structure. Note that P (a|k) implicitly depends on i and r (cf. Eq. (6)), thus being a marginal (r + 1)-body probability distribution. Again, like in the d = 2 case, by virtue of Fine's theorem [2], it suffices to consider those probability distributions in which P (a|k) factorizes; i.e., In summary, the overall minimization is performed in O(nd m(R+1) ) time. Since d, m and R are fixed, the overall scaling is O(n), although in practice it is advisable to bear in mind the prefactor. This algorithm gives not only the classical bound β C (for which it requires O(d mR ) = O(1) memory) but it also constructs a deterministic local strategy achieving it (requiring O(mn) = O(n) memory).
Let us now consider the Periodic Boundary Conditions (PBC) case. This can be reduced to the OBC case by splitting the Bell inequality with PBC at an arbitrary position i while fixing a value of M (i,R) . The correlators that contain parties in the set {i, . . . i + R} can be effectively moved to the left/right of the cut by updating the coefficients of I: If all the parties on the correlator belong to the set {i, . . . , i + R}, then this correlator has a definite value which becomes a constant in the optimization; if just one party lies outside, then the one-body weights at sites i − 1 or i + R + 1 can be updated accordingly, and so on. In Appendix B we explain this procedure in detail.
Since the amount of M (i,R) 's for which the actual minimum of I is achieved is finite, the PBC case is solved by considering d mR OBC cases. This does not change the overall complexity, but it increases the prefactor. The classical bound of I is in the PBC case found in O(nd m(2R+1) ) = O(n) time and O(1) memory for β C and O(n) memory for the deterministic local strategy.
Note that it is crucial for our method that R is constant. If R were comparable to n, the dynamic programming procedure would no longer be efficient. In fact, optimizing onedimensional Bell inequalities with full-range correlators is equivalent to optimizing general Bell inequalities, where results in computer science indicate that this is an extremely hard problem [3,24]. Note that, even in the bipartite case, deciding whether a probability distribution for m dichotomic measurements is local is NP-complete [32].

C. The quantum value
In this section, we show how to find the ground state energy of the Hamiltonian H introduced in Eq. (1). To do that, we first review the exact diagonalization of quadratic Hamiltonians in fermionic operators (Section II C 1). Then, we transform the spin operator H to a quadratic fermionic operator H via the JW transformation [33]; see also [34,35] (Section II C 2).
Throughout this paper, we will denote fermionic operators with a hat and spin operators without. Note that the JW transformation is a global operation that breaks the sense of locality in a Bell experiment, but we use it merely as a mathematical tool to find the ground state energy of H.
Recall that a quadratic Hamiltonian in fermionic operators has the form where A ij and B ij are complex numbers and h.c. stands for hermitian conjugate. Theâ i (â † i ) are annihilation (creation) operators of a fermionic system of n Dirac modes, indexed by i, with 0 ≤ i < n. The i-th mode is assigned an annihilation operatorâ i and a creation operatorâ † i . The creation operator a † i , acting on the vacuum state |Ω , populates it with a single excitation: |1 i :=â † i |Ω , whereas the annihilation operator satisfiesâ i |Ω = 0 ∀i. Such operators satisfy the following canonical anticommutation relations (CARs): The CARs (17) imply that, without loss of generality, the matrices A and B in Eq. (16) can be taken to be Hermitian and antisymmetric, respectively.

Exact diagonalization
In this section we compute the ground state energy of a Hamiltonian of the form (16). The operatorĤ of Eq. (16) can be written in terms of Majorana fermions aŝ where the 2n Majorana fermionsĉ i,α are defined aŝ where i 2 + 1 = 0. Note thatĉ i,α are Hermitian operators, and they satisfy the CARs It follows that the matrix H in Eq. (18) can be taken real antisymmetric due to Eq. (20) without loss of generality. Since every real antisymmetric matrix H admits a Williamson decomposition H = OJO T [36], where O is an orthogonal transformation and the operatorĤ is diagonalized by introducing a new set of Majorana operatorsd k,a : where {id k,0dk,1 } mutually commute and the new Majorana operators are defined aŝ achieved on a simultaneous eigenstate of the operators {id k,0dk,1 }, with respective eigenvalue s k := −sign(ε k ).

From spins to fermions
In this section we review the JW transformation and we show that the spin Hamiltonians that can be diagonalized with the method described in Section II C 1 are precisely those of the form of Eq. (1).
The JW transformation establishes an isomorphism between the Fock space of n Majorana modes and the n-qubit Hilbert space (C 2 ) ⊗n . For Majorana fermions, the JW transformation can be expressed aŝ It follows that for every fermionic operator we obtain an operator acting on (C 2 ) ⊗n and viceversa. The operatorĤ as in 2 The set of orthogonal matrices of size n is denoted O(n).
Eq. (18) consists of the terms iĉ i,0ĉi,1 , which correspond to which correspond to the string operators Str i,r α,β (cf. Eq. (2)) if i + r < n. If i + r ≥ n, we need to define a global parity operatorP In this case, the string operator takes the form Note that, since H commutes with P, its ground state has a well defined eigenvalue p = ±1 of the parity. This fact will become relevant in Section III.

III. A NEW METHOD TO OPTIMIZE BELL INEQUALITIES
The method described so far can also be seen from the opposite perspective: One considers Bell inequalities of the form of Eq. (7) with given coefficients γ (i,r) k 's and one calculates its β C with dynamic programming. On the other hand, we optimize the quantum value over a restricted set of measurements, namely, single qubit σ z measurements, and for the r-body correlators, qubit measurements in the x − y plane for parties i and i + r and σ z measurements in the intermediate ones (cf. Eq. (7)). The resulting Bell operator can be mapped to a system of fermions as in Eq. (18), whose diagonalization allows us to find the minimal quantum value.
The latter process is carried out as follows. Let us denote the k-th observable of the i-th party by M (i) k , with k ranging from 0 to m − 1 (or m if R > 1). We shall pick qubit observables parametrized as M y (k < m) and construct the Bell operator B, which will be of the form (4). Note that to build the Bell operator, one simply needs to substitute the correlators in (7) by the corresponding quantum observables: If there exists a quantum state ρ for which Tr(Bρ) < 0, then ρ is nonlocal. If this is the case, we shall denote the quantum violation observed by QV := Tr(Bρ). Note that ρ can be taken, without loss of generality, as a projector onto a ground state of H. We will look for the optimal measurement settings ϕ (i) * k such that B has the minimal eigenvalue. The spin Hamiltonian H is again diagonalized by applying the JW transformation (cf. Eq. (25)), and Eq. (1) is almost transformed into Eq. (18) up to the string operators that cross the origin, which carry a parity operatorP. Since the eigenstate with the lowest eigenvalue ofĤ has a well defined parity p = ±1 (because [Ĥ,P] = 0), we can change P by p in Eq.
(28) so thatĤ is now quadratic. One has to make sure, though, that the superselection rule imposed by initially choosing p is obeyed. That is, the ground state ofĤ needs to satisfy Eq. (30) stems from the fact that, under the transformations of Eq. (23),P is transformed as (see Appendix A for a proof) If Eq. (30) does not hold, one has to modify Eq. (24) accordingly by picking The minimal E 0 for p = 1 or p = −1 is the minimal eigenvalue of B.
Finally, note that if R = 1 and m = 2, the minimal eigenvalue of B(ϕ (i) * k ) yields the minimal value of I achievable within quantum theory, denoted −β Q , because the optimal quantum violation of Bell inequalities with n parties and two dichotomic observables is obtained with qubits and traceless observables on a plane [37]. However, if R > 1 or m > 2 this result does not hold in general, as higher-dimensional systems and more general observables can produce more nonlocal correlations.

IV. THE TRANSLATIONALLY INVARIANT CASE
In this section, we consider the case in which H in Eq. (1) is translationally invariant (TI). In the spirit of Section III, the following procedure can be seen as an optimization of a TI Bell inequality with the same set of observables at each site. We first present an algorithm to compute the classical bound of a TI Bell inequality with short range correlators which is exponentially faster in the number of parties than the one of Section II B 1 (Section IV A). We then find the ground state energy of a TI Hamiltonian of the form of Eq. (1) analytically (Section IV B).

A. Exponentially faster solution of the classical bound
In this section, we present a method that is exponentially faster in the number of particles than the one in Section II B 1 to compute the classical bound of TI Bell inequalities of the form of Eq. (8); i.e., where γ (i,r) k are independent of i. For the sake of simplicity, we first present this problem in a more abstract setting. Then we adapt it to TI Bell inequalities.
Consider a function f (0) : S × S −→ R where S is a finite set. We will describe how to compute in O(log 2 w) steps. To this end, let us define for t > 0. Note that the superscript t indicates the iteration step. The idea is to successively rewrite F in terms of f (t+1) instead of f (t) thus eliminating at each step approximately half of the variables in the minimization. For instance, by writing F in terms of f (1) one has already carried out the minimization over all the variables with odd index (except the last one if w is odd). Note that any function f (t) is defined by specifying |S| 2 numbers, where |S| is the number of elements of S. Computing Eq. (34) requires O(|S| 3 ) operations. Note that when w is a power of 2, then F is given by In the general case, however, we cannot assume w to be a power of 2. Nevertheless, every positive integer w can be uniquely expressed as a sum of different powers of 2. The idea is to apply the procedure described above to each of these powers of 2 and then optimize over the remaining O(log 2 w) variables. To this end, recall that the binary expression for w is where · is the floor function, a i ∈ {0, 1} correspond to the digits of w in binary, |w| := i a i is the Hamming weight of w, and the b j 's enumerate the indices i for which a i = 1, sorted in decreasing order. For instance, if w = 11, (a 3 , a 2 , a 1 , a 0 ) = (1, 0, 1, 1) and (b 0 , b 1 , b 2 ) = (3, 1, 0), and if w = 15, then (a 3 , a 2 , a 1 , a 0 ) = (1, 1, 1, 1) and Note that F can now be rewritten as which is a minimization over 1 + |w| ≤ 1 + log 2 (w) variables y j , where · is the ceiling function. Note that the f (bj ) need no longer be the same function for different j, so the expression for F given in Eq. (37) is not TI in general. Note that at the log 2 w -th step of the optimization we have the x 0 variable on the left (see Fig. 2) and several remaining variables x j for j ≥ log 2 w . Since these remaining variables highly depend on the binary expression of w, we denote them by y j (cf. Eq. (37)). To minimize over y j , we proceed in a similar fashion, now defining g (0) := f (b0) and x 0 Here we find F for w = 11. We represent the variables xi with circles and the functions f or g with lines. Each line with the same color corresponds to the same function. The variables that are not in gray are eliminated at the next step. At the 0-th iteration, there are 10 functions f (0) depending on 11 variables. We compute f (1) and we substitute it as many times as possible, so that at the 1-st iteration, we have eliminated the variables in purple. Then we eliminate the variables in blue by computing f (2) and the variables in turqouise by computing f (3) . We then compute the functions g by joining them with the remaining f 's, thus eliminating the orange and the red variables. Finally we minimize on the ends, where we can impose conditions on the boundary if needed.
for t > 0. Observe that now the optimization is linear, similar to the dynamic programming presented in Section II B 1. It follows that F can now be computed as In Figure 2 we describe this procedure with an example.
To adapt this algorithm to the minimization of I, we start by noting that I can be written as where h is defined as and the indices of the parties are taken modulo n in the M (i,r) k defined in Eq. (9). Observe that every i-th and (i + R + 1)-th parties share R parties via h. In particular, by picking i = j(R + 1) − 1 we denote their local deterministic strategy as We now rewrite the optimization of I over M in terms of x j . To this end, let us define x 0 3. An example with n = 14 and R = 3 (this corresponds to w = 3). Each circle or square corresponds to a party, starting at i = 0 on the left. The dash-dotted line represents crossing the origin and each line below the parties represents a function h (note that h has a range of R + 1 parties). The lines corresponding to the h's are arranged in groups of R + 1 (except for the T corresponding to the tail), which we represent with the same color. Each full group contains a single gray square. By minimizing the local deterministic strategy at the square, we can define f (0) which depends on the local deterministic strategy chosen at the R neighbours of each side. These groups of R parties correspond to the xj in Fig. 2. They encode the possible d mR local deterministic strategies for each xj. Now, in O(log n) steps we find g |w|−1 (x0, xw). Finally, we find the classical bound by minimizing the sum of g |w|−1 and T restricted to the x0 and xw that have a compatible overlap.
Since we cannot assume that n is a multiple of R + 1, we take w := n/(R + 1) . Then, where the tail T (x w , x 0 ) is defined as with the indices of the parties taken modulo n and min is the minimum taken on those x 0 , x w that are compatible with PBC (for instance, if w = n(R + 1), then T = 0 and min is taken over x 0 = x w ). In Fig. 3 we illustrate the procedure we described with an example.

B. Analytical solution of the quantum value
Here we consider Eq. (1) in the TI case, which corresponds to t (i) and t (i,r) α,β being independent of i. In terms of Bell inequalities, this corresponds to the optimization of TI Bell inequalities with the same set observables being performed at each site. We give analytically closed expressions in this case. As we prove in Appendix C, the Williamson eigenvalues for a TI Bell operator of the form of Eq. (4) are given by with k ranging from 1 to (n − (p − 1)/2)/2 , where with υ k,r := rπ(2k − (p + 1)/2)/n. Depending on the parity of n, the following eigenvalues also appear: If p = −1, then ε 0,+ always appears and ε 0,− only appears if n is even. If p = 1, then ε 0,− only appears if n is odd and ε 0,+ does not appear (see Appendix C). The superselection rule (30) to be fulfilled in this case reads Note that if we just want to find the ground state energy of a TI fermionic Hamiltonian (16), then the matrices A and B are circulant, which means that the previous analysis can be done assuming that p = −1 and no superselection rule needs to be obeyed in this case, as Eq. (52) appears from the transformation of spins to fermions. This result is applied to the example of Section V A.

V. EXAMPLES
In this section we present three different examples in which we illustrate the tools we have presented through the paper. In Example V A we optimize a tight TI inequality for n = 8 parties with R = 2 and PBC, showing that it has quantum violation when the same set of qubit measurements are performed at each site. Interestingly, such optimal measurements are M 0 = σ x , M 1 = σ y and M 2 = σ z . In Example V B we construct a quasi TI Bell inequality for any even number of parties and any number of measurements, which depends on one parameter. We find its classical bound analytically with dynamic programming. The Bell operator corresponds to a spin XY model which we also solve analytically. Finally, in Example V C we show that the ground state of a spin glass is nonlocal in some parameter region.

A. A translationally invariant Bell inequality for 8 parties
The general form of a translationally invariant Bell inequality I + β C ≥ 0 with m = d = R = 2 is (cf. Eq. (7)) where the translationally invariant correlators T are defined as In Appendix D we present a table with the optimal (tight) Bell inequalities of these kind for n ≤ 8 and the quantum violation we can observe. Let us remark that finding all Bell inequalities for a given scenario is a computationally very expensive task and one typically manages to do it only for very small values of n, m, d and R, even if symmetries are imposed [38]. When looking for Bell inequalities of the form (53), n = 8 was the maximum number of parties for which this task could be carried out in a reasonable time. Here we present a particular case as an example. If one takes the following coefficients: γ = 0, γ 00 = γ 10 = −γ 01 = −γ 11 = 2, −γ 020 = −γ 021 = γ 120 = γ 121 = 1, then the dynamic programming gives β C = 32 and the mea- The latter is proven by applying Eq. (46) to the example. More specifically, we observe that the chosen coefficients and measurements yield an H matrix (cf. Eq. (18) and Appendix C) with upper-diagonal blocks h 1 and h 2 and the rest of the h r blocks are zero. This greatly simplifies the expression for ε k as x k = 0 and a k = c k imply ε k = 2(a k ± |b k |) (Note, that in the range of interest, b k ≤ 0 so that ε k = 2(a k ∓ b k )). If we take the plus sign, we have and if we take the minus sign we obtain ε k,− = 16 cos π 4k + 11 − p 16 .
Thus, we can now calculate the ground state energy consistent with each p, which is given by if p = −1, and E 0 = −8 √ 2 + 2 cos(π/8) + 2 sin(π/8) if p = 1. One does not need to check the superselection rule Eq. (52) for p = −1 as some of the ε k are zero. However, we do need to check it for p = 1. There, we took an even (2) number of sign flips to the ε k and the determinant of O is 1 (cf. Appendix C) . It follows that Eq. (30) holds. The case for p = −1 gives E 0 −27.3137 whereas for p = 1 it gives E 0 = −32.2187. Hence, I + β C = −0.2187 < 0, signalling the presence of nonlocality.

B. A quasi translationally invariant Bell inequality
Let us consider the chained Bell inequality [39] between two parties labelled A and B: where I (A,B) chain is given by where it is assumed that A −1 := −A m−1 . Note that the CHSH inequality [40] I

is a particular case of Eq. (59) for m = 2. Inequality (58) is violated maximally with the following settings
and B i = A i , where the angles are given by φ i := (i+1)π/m, and with the state giving β Q = I (A,B) chain = −2m cos(π/2m). Notice that the maximal violation relative to the classical bound is β r Q := β Q /β C = [m/(m − 1)] cos(π/2m).
The bipartite Bell operator corresponding to the chained Bell inequality with the above measurements can be written as where α m := m cos 2 (π/2m) and β m := (m/2) sin(π/m). By defining σ π/2m := cos(π/2m)σ x +sin(π/2m)σ y , this operator can be further re-expressed in a formally similar manner to the XY Hamiltonian as Let us now consider the following Hamiltonian: where the weights f i ( ) alternate from even to odd sites as f i ( ) := 1 + (−1) i and is an arbitrary real parameter. We note that the Hamiltonian (64) is a particular case of the onedimensional Bell inequality when the same measurements (60) are taken at each site.
Let us now determine the classical bound of I (2n) chain ( ). For even 3 n, the dynamic programming gives β C = 4n(m − 1) max{1, | |}. The explanation for this result is that, whenever ε > 1, it is better to use a classical strategy that would give −2(m − 1) on every I (2i,2i+1) chain inequality and 2(m − 1) on every I (2i+1,2(i+1)) chain inequality. An exemplary local strategy achieving this bound is given by It is worth highlighting two limiting cases: • = 1. This corresponds to a sum of disjoint chained Bell inequalities between pairs (2i, 2i + 1).
which is maximally violated by the state |ψ m AB ⊗ |ψ m CD ⊗ · · · . The quantum value is then β Q = 2nm cos(π/2m). Hence, there is a O(1) violation relative to the classical bound that holds for every n: • = 0. This case corresponds to a sum of the chained Bell inequalities with the same weights between neighbours: This inequality cannot be violated, as quantum correlations are monogamous with respect to the chained Bell inequality [37,41]. Loosely speaking, if party B violates I chain with A, it cannot violate it simultaneously with C. For some types of monogamy relations, this result holds for various generalizations of the CHSH inequality to more measurements, outcomes and parties [42].
It is then clear that there is some critical value of for which correlations stop being nonlocal and one is able to simulate them locally.
Let us now notice that the 4n × 4n matrix H appearing in Eq. (18) and corresponding to the Hamiltonian (64) has the form H = H 0 ⊗ H 1 , where with f 0 := 1 + and f 1 = 1 − for short, and This is seen by applying the JW transformation to Eq. (62). In this case, to find the Williamson eigenvalues of H it is sufficient to find those of H 0 , which will appear with both signs, as H 1 has eigenvalues ± cos(π/2m).
Hence, the quantum bound will be β Q ( ) = 2 cos π 2m Although one should check the superselection rule Eq. (30), it is not necessary for large n, as the difference between β Q for p = 1 and β Q for p = −1 vanishes as n grows. Let us analyze the behaviour of β Q ( ) in the thermodynamic limit. The contribution per particle to β Q ( ), denoted which is a Riemann sum, so it is by definition This can be expressed more compactly as where E(t) is the complete elliptic integral of the second kind 4 . In Figure 4 we can see such behavior compared to the contribution per particle to the classical bound,β C ( ) = 2(m − 1) max{1, | |}. 4 The complete Elliptic integral of the second kind is defined as with the parameter t obeying 0 < t < 1. The plot corresponds to a spin glass of n = 100 spins with PBC, averaged over 1000 realizations. The horizontal axis corresponds to the mean µ of the Gaussian distribution and the vertical axis corresponds to its standard deviation σ. If µ = 0, the value is constant for all σ > 0, as both the expected value of the ground state and the classical bound grow at the same rate with σ. The white line represents the level curve for βC = βQ. The top-left part of the plot corresponds to the region of parameters in which we detect nonlocality (βC < βQ). Note that on the bottom-right region one finds values for which βQ/βC < 1 due to the fact that there is no simultaneous eigenvalue of σ π/4 and σy (the classical bound cannot be saturated using σ π/4 and σy as observables).

C. A spin glass
Finally we consider a Hamiltonian similar to Eq. (64) for m = 2, where all the couplings are random, generated from a Gaussian probability distribution with mean µ and standard deviation σ, like a spin glass: We can efficiently compute the ground state of Eq. (75) and the classical bound of the CHSH-like Bell inequality associated to it with the methods we have presented, although numerically. We expect that, when |µ|/σ 1, we do not detect nonlocality, as the resulting inequality is close to the monogamous limit = 0 of Example V B. However, if σ is big enough, one expects to have links with high weight surrounded with links with smaller weight; then, it may compensate to violate those links with higher weight. In Figure 6 we show the result, for n = 10 2 spins, averaged over 10 3 realizations. There is clearly a transition for which the complexity of the ground state allows or does not allow for the statistics that one obtains when performing measurements on it to be simulable locally.
D. An XXZ-like spin model based on Gisin's Elegant Bell inequality In this last example, we present an XXZ-like Hamiltonian, which is not solvable via the JW transformation. We find its ground state energy numerically, using tensor networks and DMRG [30]. In this case we find a much richer structure than in Example V B. The Bell inequality that we associate to this Hamiltonian is a modification of Gisin's Elegant Bell inequality [29]. Gisin's original inequality is defined in a bipartite scenario with four dichotomic measurements with outcomes ±1 on Alice and three dichotomic measurements with outcomes ±1 on Bob: and it reads |I| ≤ 6. We observe that with the following observables: the corresponding Bell operator becomes and its ground state is |ψ − = (|01 − |10 )/ √ 2, yielding an expectation value ψ − |B|ψ − = −4 √ 3 −6.9282. Let us now introduce the following modification, where ∆ is a real parameter, where S ∆ is a 4 × 3 matrix defined as The classical bound becomes Now we have that the Bell operator has become in either case. Its ground state energy is The ground state is |ψ − if ∆ > −1 and it lies in the subspace spanned by |00 and |11 if ∆ ≤ −1.
We can now construct the many-body Bell inequality in a similar fashion as in Example V B: Note that the Bell inequality J has 4 binary measurements with outcomes ±1 on the even sites and 3 binary measurements with outcomes ±1 on the odd sites.
The dynamic programming procedure yields the following classical bound in terms of and ∆, which is a C 0 piece-wise continous function. Due to all the cases that appear, and the complexity of the inequality, we omit the description of the local deterministic strategy. The regions are defined as follows (see Fig. 7): If n ≡ 2 mod 4, n > 2, the classical bound is, on each region: whereas if n ≡ 0 mod 4, the classical bound simplifies to Using the A k measurements on the even sites and the B l ones on the odd sites, this yields the following XXZ-type Hamiltonian: The ground state of Eq. (89) does not have an analytically closed form and has to be computed using DMRG. To do so, we have used the iTensor library [30]. The results are plotted in Fig. 8, where we observe quantum violations in a parameter region that does not seem to vanish as n grows. We also observe a different behavior depending on the parity of n/2, which we attribute to some sort of frustration arising in the classical optimization, especially for low values of n.

VI. CONCLUSIONS AND OUTLOOK
In this work we have shown, on the one hand, that the ground states of some spin Hamiltonians are nonlocal. We have associated a Bell inequality to the Hamiltonian and we have computed its classical bound. We have exactly diagonalized the Hamiltonian and we have computed the quantum value of the corresponding Bell inequality. We have achieved this goal by combining two rather unexplored techniques in  Fig. 7, also sketched here in gray (for odd n/2) and black (even n/2) for clarity. The region for which there is quantum violation is the bounded set (i.e., containing the point (∆, ) = (1, 1)). The plot is symmetric with respect to the line = 0.
quantum information, which are dynamic programming and the Jordan-Wigner transformation. This has allowed us to detect the presence of nonlocal correlations. On the other hand, these tools have provided us a new way to determine the classical bound of certain classes of Bell inequalities, and look for their quantum violation under conveniently chosen observables. In the case of two dichotomic measurements, the optimization also yields the quantum bound (the maximal quantum violation) of the inequality. In the TI case, we have provided an algorithm to find the classical bound that is exponentially faster in the system size and we have obtained exact analytical solutions for the quantum value. Then, we have applied these techniques to several examples in which we reveal nonlocality: A tight TI Bell inequality with a TI Bell operator for 8 parties, a quasi TI Bell inequality with a uniparametric quasi TI XY Hamiltonian and the ground state of a XY spin glass in some parameter region. We have also seen that for interacting models such as XXZ-like spin Hamiltonians, our method can be applied. There, we find the ground state energy numerically and we map the Hamiltonian to a modified version of Gisin's Elegant Bell inequality. These findings open the possibility to the implementation of multipartite quantum information protocols that use nonlocality as a key resource, by using ground states of Hamiltonians that appear naturally. We remark that the Hamiltonians and the Bell inequalities we have studied have a finite interaction range, which in the TI case makes them particularly interesting from an experimentally-friendly perspective. Note that previous Bell inequalities for quantum many-body systems with low order correlators were specially designed for a permutationally invariant (PI) symmetry [25,26]; while being able to detect nonlocality in ground states of physical Hamiltonians such as the Lipkin-Meshkov-Glick [43] or a spin-squeezed Bose-Einstein condensate [27], the information accessible to these inequalities is bound to a de Finetti theorem [44,45], thus becoming more compatible with that produced by a separable state as the system grows [26]. This requires to increase the number of measurements with the system size in order to close the finite-statistics loophole [27]. However, many systems of interest are not PI, but TI, and we have studied the latter in this work. In this case, a de Finetti restriction does not apply, making the detection of nonlocality more robust to experimental imperfections. Interestingly, to study the Bell inequalities proposed in [25,26], in the classical and quantum many-body regime, one employs powerful mathematical tools (namely, convex hulls of semialgebraic sets [46] or the Schur-Weyl duality from representation theory [47], respectively), which no longer apply when the PI symmetry is broken. Here we have used other mathematical tools (dynamic programming for the classical bound and the JW transformation for the quantum value) that allowed us to solve, even exactly, the two cases with this much weaker symmetry.
Let us finalize by pointing out possible future research directions that stem from our work. Throughout this paper, we have restricted ourselves to the study of nonlocality in one-dimensional spin short-ranged Hamiltonians, as we could compute their ground state energy with exact diagonalization. One can eliminate this restriction and study short-ranged spin Hamiltonians by using a Matrix Product State (MPS) description of the state, which is a good ansatz for these systems [48]. The Bell inequalities that we would naturally associate to them would still be of the form of Eq. (8), because the interaction range R would be fixed, so we could still efficiently find their classical bound. Conversely, we can eliminate the restriction on the subset of observables that we choose in Section III and the string of σ z 's in the middle of the string operators, thus studying purely two-body correlator inequalities like the classes derived in [38]. In such cases, powerful numerical algorithms such as DMRG [49] would be suitably tailored to perform the quantum optimization and check whether nonlocality is detected. Another interesting problem is related to the persistency of nonlocality [50]. While the tools to carry this kind of analysis in the PI case have been put forward [26], in the case of one spatial dimension we do not know yet how robust are the inequalities we have presented to particle losses. One could also generalize our results towards other directions. Since Bell inequalities with more than two inputs or outputs per party are not, in general, maximally violated by measuring qubits [51], if one would like to increase the chance to detect nonlocality with these more general classes of Bell inequalities, increasing the physical dimension of the system would be a way to obtain better results. The tools presented here could be extended by using a generalized JW transformation [52] from qudits to parafermions, although the problem becomes algebraically more involved.
We have also seen through the examples presented that there is a strong relation between Hamiltonians of physical systems and Bell inequalities that we associate to them. Whereas we naturally establish this connection in the following direction: starting from the Hamiltonian of a physical system, we assign a Bell inequality to it, one can think of this relation in a more general way, since there are many Bell inequalities that, with the appropriate observables, realize the same physical Hamiltonian. Moreover, in Example V D, we have seen that this correspondence can be non-trivial, as the inequality does not even have the same number of measurements at every site. With the tools we have presented here, it is now possible that, given a physical one-dimensional Hamiltonian with short-range interactions, tailor the best Bell inequality that reveals the nonlocality of its ground state, as such Hamiltonian corresponds to the realization of a Bell operator of many different Bell inequalities, each with its own classical bound, which we can compute efficiently. Proof. To prove Eq. (A5), we can split the sum into the indices for which (i, α) = (j, β) and the indices for which (i, α) = (j, β).
In the second case, we note that we can rewrite the sum as u k,0 u k,1 because of the CARs (20).
In order to prove Eq. (A6), we begin by noting that We can now split the sum (A6) into the parts where i = k and i = k.
For the first part, we have a contribution in Eq. (A6) that amounts to i n−1 k=0 (u 2 k,0 + u 2 k,1 )ĉ k,0ĉk,1 + 2u k,0 u k,1 1 |Ω , which, using Eq. (A8), simplifies to For the second part, the contribution to Eq. (A6) is expanding the sum over α we have i =k Splitting the sum between those indices for which i < k and those for which i > k we have that it can be rewritten into an expression involving the CARs (17): where the product is written from left to right. Note that n k can only be 0 or 1, due to (17). The CARs (17) further imply thatâ i |N = (−1) i−1 j=0 nj |N whenever n i = 1 (N and N differ only in its i-th index) andâ i |N = 0 whenever n i = 0.
Let |ψ be an eigenvector ofP c of eigenvalue 1. We can assume, without loss of generality, that |ψ is a basis element |N of the Fock space, with n−1 i=0 n i ≡ 0 mod 2. By noting that id k,0dk,1 can be expressed as we see that the only operators which appear are products of a i orâ † i withâ j orâ † j . Hence, they either annihilate |N or they add −2, 0 or 2 to its particle number; always conserving its parity. This leaves us with a linear combination of vectors that live in the same subspace F e , so it remains on F e . The same argument applies to F o .
Hence, the subspaces F e and F o of the operatorP c are invariant under orthogonal transformations of the Majorana fermions. However, underP d , they might change its parity. Because of linearity and Eq. (A13), it suffices to prove that |Ω is an eigenstate ofP d , as it stems from Lemma 2. If its eigenvalue is 1, then F e = F e ; if it is −1, then F e = F o .
Appendix B: The PBC to OBC reduction In this appendix we describe the reduction of the optimization problem of finding the classical bound in a Bell inequality with PBC to the optimization problem of finding the classical bound for a Bell inequality with OBC in one dimension.
For the sake of clarity, we first describe the procedure for inequalities with an arbitrary interaction range R, but without the chain of M (j) m observables in the middle of the string operators (cf. Eq. (7)), so that the inequalities strictly contain one or two-body correlators. We pick R consecutive parties. Without loss of generality, we can assume them to be labelled from 0 to R − 1. To these parties we assign one of the d m·R possible deterministic local strategies. Let (i, j) be a pair of parties. There are three cases to consider: • If i and j are between 0 and R − 1, any correlator between parties i and j has now a definite value.
• If either i or j, but not both are between 0 and R − 1, only one side of the correlator has a definite value. Because the classical bound is computed on a deterministic local strategy, the value of the correlator factorizes as the product of the local values assigned by the strategy we are considering. Hence, these correlators can be effectively moved outside of the interval 0 . . . R − 1 by updating the one-body term of the party outside that interval.
• In any other case, the correlator remains the same.
In Figure 9 a) and Figure 9 b) we illustrate this procedure for an example with R = 2. The first case is illustrated in blue; the second is illustrated in red and the last one is illustrated in black.
Finally, if we have an inequality of the form of Eq. (7), the intermediate M m values closest to the boundary become fixed (2(R−2) of them). This has to be taken into account when performing the dynamical programming, as one should only explore those configurations compatible with the boundary conditions imposed by the PBC problem with the deterministic local strategy that we have chosen to make the reduction.

Appendix C: Quantum optimization of translationally invariant Bell inequalities
In this appendix we discuss the quantum optimization for the translationally invariant (TI) case (i.e., when the Bell in-Real DFT (RDFT) matrix of order n: (R n ) kl := 2 n cos 2πkl n − π 4 , 0 ≤ k, l < n. (C3) It is easy to see that R n = R T n and R 2 n = 1, so R n is orthogonal. In the following, we are going to use the fact that for any n, det(R n ⊗ 1 2 ) = (det R n ) 2 = 1.
Let us now, with the aid of the orthogonal transformation R n , study which invariant subspaces H acts upon. If p = −1, then we compute (R n ⊗ 1 2 )H(R n ⊗ 1 2 ); if p = 1 then we calculate (R 2n ⊗ 1 2 )H(R 2n ⊗ 1 2 ). Note that 1 2 acts on the two Majorana modes asociated to one site. On which corresponds to the Williamson eigenvalue(s) for the 2× 2 blocks.
To bring H to a Williamson form, it remains to bring the 4 × 4 blocks J k and J k to a Williamson form. To this end, let us start by defining υ k,q := qπ(2k − (p + 1)/2)/n and Let us notice that the blocks J k (or J k ) take the following form:  Now we are ready to find an orthogonal transformation that brings J k or J k to a Williamson form, which we state in the following lemma: Its two Williamson eigenvalues ε k,± , are given by where ∆ k := (a k − c k ) 2 + 4(b 2 k + x 2 k ) and k ranges from 1 to n/2 + (p − 1)/4 . Furthermore, this orthogonal transformation always satisfies det O k = −1.
Proof. The choice of O k is not unique in general. Here we are going to construct O k as the product of three matrices. Since k is fixed, we are going to skip explicitly stating the subindex throughout the proof. We construct O as O := LM R, where L and R are diagonal matrices whose entries are defined by: and R −1 = 2 b 2 + x 2 Diag{r − , r − , r + , r + }, where r ± := ∆ ± (a − c) √ ∆. The matrix M in the middle is not diagonal, and in its most general form, can depend on two real parameters, which we denote φ and θ. We have that the entries of M are given by  Tracking all the transformations we have made, i.e., counting the parity flips imposed by the choice of all the orthogonal transformations we have made, the superselection rule Eq. (30) that must be obeyed takes the form of Eq. (52).
If, in addition, we impose a finite interaction range R, then we can further simplify the expressions for x k , a k , b k and c k thanks to the property h r = −h T n−r and arrive at Eqs. (47,48,49,50).