Unveiling quantum entanglement in many-body systems from partial information

Quantum entanglement is commonly assumed to be a central resource for quantum computing and quantum simulation. Nonetheless, the capability to detect it in many-body systems is severely limited by the absence of sufficiently scalable and flexible certification tools. This issue is particularly critical in situations where the structure of entanglement is a priori unknown, and where one cannot rely on existing entanglement witnesses. Here, we implement a scheme in which the knowledge of the mean value of arbitrary observables can be used to probe multipartite entanglement in a scalable, certified and systematic manner. Specifically, we rely on positive semidefinite conditions, independent of partial-transposition-based criteria, necessarily obeyed if the data can be reproduced by a separable state. The violation of any of these conditions yields a specific entanglement witness, tailored to the data of interest, revealing the salient features of the data which are impossible to reproduce without entanglement. We validate this approach by probing theoretical many-body states of several hundreds of qubits relevant to existing experiments: a single-particle quench in a one-dimensional $XX$ chain; a many-body quench in a two-dimensional $XX$ model with $1/r^3$ interactions; and thermal equilibrium states of Heisenberg and transverse-field Ising chains. In all cases, these investigations have lead us to discover new entanglement witnesses, some of which could be characterized analytically, generalizing existing results in the literature. In summary, our paper introduces a flexible data-driven entanglement detection technique for uncharacterized quantum many-body states, of immediate relevance to experiments in a quantum advantage regime.

Quantum entanglement is commonly assumed to be a central resource for quantum computing and quantum simulation. Nonetheless, the capability to detect it in many-body systems is severely limited by the absence of sufficiently scalable and flexible certification tools. This issue is particularly critical in situations where the structure of entanglement is a priori unknown, and where one cannot rely on existing entanglement witnesses. Here, we implement a scheme in which the knowledge of the mean value of arbitrary observables can be used to probe multipartite entanglement in a scalable, certified and systematic manner. Specifically, we rely on positive semidefinite conditions, independent of partial-transposition-based criteria, necessarily obeyed if the data can be reproduced by a separable state. The violation of any of these conditions yields a specific entanglement witness, tailored to the data of interest, revealing the salient features of the data which are impossible to reproduce without entanglement. We validate this approach by probing theoretical many-body states of several hundreds of qubits relevant to existing experiments: a single-particle quench in a one-dimensional XX chain; a many-body quench in a two-dimensional XX model with 1/r 3 interactions; and thermal equilibrium states of Heisenberg and transverse-field Ising chains. In all cases, these investigations have lead us to discover new entanglement witnesses, some of which could be characterized analytically, generalizing existing results in the literature. In summary, our paper introduces a flexible data-driven entanglement detection technique for uncharacterized quantum many-body states, of immediate relevance to experiments in a quantum advantage regime.

I. INTRODUCTION
Quantum entanglement is a distinguished feature of composite quantum systems, marking a fundamental departure from their classical counterparts [1]. Over the last decade, it has become a commonplace that manybody entanglement represents an essential resource for quantum computation [2], quantum simulation [3], and quantum metrology [4]. While, on the theoretical side, the exact role of quantum entanglement in offering a quantum advantage remains somewhat controversial [77], the ability to manipulate quantum many-body superpositions arguably represents a major endeavour for many experimental platforms. As a matter of fact, the controlled preparation of many-body entangled states is a hallmark of such capability, and has been achieved in several experimental systems [5][6][7][8][9][10][11][12][13]. On the other hand, a growing number of experiments operate in regimes inaccessible to the best available classical simulations [14][15][16][17][18][19][20]-another hallmark pointing towards a genuine quantum advantage. It is commonly assumed that the intractability of classical simulations originates in the largescale quantum entanglement which develops across the experimental system [2]. Nevertheless, a proper quantum computation, performed in a regime inaccessible to the best classical algorithms, and where the structure of quantum entanglement is also probed, has not yet been reported. This absence is partly due to the lack of sufficiently flexible and scalable theoretical tools to analyze the experimental data produced in such quantum devices. Surely, one cannot simply rely on the violation of existing entanglement witnesses, for the structure of entanglement in the system, and therefore the suitable entanglement criterion to potentially reveal it, are a priori unknown. Furthermore, one cannot envision to use tomographic information about the underlying quantum state -for acquiring such information would require a number of measurements growing exponentially with the system size [21,22].
It is precisely the purpose of the present paper to show the broad applicability of a flexible and scalable tool to certify entanglement in an unknown multipartite quantum state. In our setting, we assume that the expectation values of (a scalable number of) arbitrary observables are known. Starting from the same insight as in [23], namely the connection between the compatibility of data with (quantum) separable states and the (classical) truncated moment problem, we provide a simple and scalable method to construct an entanglement witness from the observed expectation values that is tailored to be robust to noise. The main insight from [23] is that if the underlying quantum state is separable, then the available data are obtained as entries of a correlation matrix, which satisfies certain positive-semidefinite constraints. Such compatibility conditions can be efficiently verified via socalled semidefinite-programming (SDP) techniques [24], allowing the study of systems of hundreds of qubits. The arXiv:2107.03944v3 [quant-ph] 16 Mar 2022 failure for the data to pass this test serves directly as an entanglement detection method, in which case our approach delivers a specific entanglement witness, violated by the observed data.
The resulting method is platform-agnostic, in the sense that how such data should be a priori chosen, and how they should be inferred in an actual experiment is not relevant, and, in fact, will not be discussed in this work. As an illustration, we consider one-and two-body correlations for N qubits, but our scheme is flexible and can incorporate the knowledge of any k-point function, or in general of any many-body observable. By benchmarking the method on paradigmatic quench experiments, we show its wide applicability and its ability to extract physically relevant entanglement witnesses. The expression of the witnesses themselves provides qualitative insight into the driving mechanism responsible for entanglement within the system. As a matter of fact, for several of the examples we have considered, we could analytically characterize the witnesses obtained numerically. This led us to extend some known entanglement criteria of the literature and derive completely new ones as well. Analagously to the hierarchy introduced in [23], the scheme we propose can be generalised as a complete hierarchy of positive-semidefinite tests: if no separable state can reproduce the available data, the data will necessarily fail to pass the test at a finite level of the hierarchyin this sense, the hierarchy is complete.
Comparison to previous works. A large body of literature has already considered the problem of entanglement detection from partial information. Some of these results are recovered as special cases of the approach implemented in this paper; some others lack the scalability required to apply them to many-body systems; and some altenative scalable schemes either lack the flexibility of the present approach, or can be inconclusive. In particular, the so-called covariance matrix criterion [25] and the generalized spin-squeezing inequalities [26], which are based on one-and two-body correlations, are recovered as a consequence of our approach (as further discussed in Appendix C) -while our approach is more flexible, as it can incorporate the knowledge of any correlation function. Criteria based on higher-order correlations have also been derived [27][28][29]. However, the efficiency of these approaches is unclear if only partial information is available (for example, if only two-body correlations are known). Furthermore, these approaches [25][26][27][28][29] provide only sufficient conditions for entanglement, and therefore can be inconclusive even though the available data cannot be reproduced by a separable state; in contrast, here we provide a systematic and convergent hierarchy of criteria. A systematic approach, which can also incorporate partial knowledge about the quantum state, was proposed based on the solution of so-called separability eigenvalue equations [30,31]; but this approach has an exponential cost and cannot be applied already to a few tens of qubits. Conceptually-different approaches, based on randomized measurements, have also been developed.
Such approaches allow one to test bipartite entanglement criteria based on Rényi entropies [12,22], and PT-based criteria [32,33]. However, in addition to the very high experimental requirements underlying these approaches, they require a number of measurements scaling exponentially with N , severely limiting their scalability beyond a few tens of qubits. Recently, intrinsically scalable approaches to the problem of multipartite entanglement detection from partial information have been developed. An entanglement-detection method from the knowledge of two-body reduced density matrices was developed in Ref. [34], with a similar computational cost as the one in the present work; however, the above approach lacks the flexibility to be adapted to an arbitrary set of data, especially the average value of many-body observables. The approach of the present paper is complementary to Ref. [35], where the problem is solved through a mapping onto an inverse problem of classical statistical physics, offering a systematic and scalable solution; however this approach could be inconclusive for particular data sets. In contrast, here we solve a relaxation to this problem with an efficiency which is independent of the structure of the data, obtaining entanglement witnesses whose violation is guaranteed by semidefinite-positive constraints. Lastly, our approach shares some similarities with the method presented in [36]. However, in contrast to ours, the method in [36] is device-independent, namely, it exploits no information about the underlying Hilbert space. This makes the resulting entanglement test sensitive to a careful choice of measurement basis for each particle.
The first entanglement detection approach based on the connection to the classical moment problem was introduced in Ref. [23]. What we develop here can be seen as a complementary separability test, based on the same conceptual premises, by with a different physical motivation. While [23] aims at constructing -if it existsa separable state compatible with the data, our focus is instead on building a scalable and robust criterion for entanglement based on the available data. Technically, we define the SDP as a noise robustness problem, and restrict our analysis to the first level of a hierarchy of conditions in order to preserve scalability. In contrast, the hierarchy in [23] exploits a cost function tailored to identify a so-called "flat extension", which is a property that, by definition, can be assessed only by solving SDPs of increasing levels in the hierarchy. However, it should be emphasized that increasing by just one level the hierarchy is already extremely costly in a many-body setting, and it is not doable already for systems of few tens of particles. Moreover, while an entanglement witness could potentially be obtained from the dual of the first level of the hierarchy in [23], it has no guarantee to be robust against noise. Lastly, we notice that the element of randomness in the objective function in [23] implies that the witness will be different for every run of the SDP, while our benchmarks allows one to derive analytical witnesses in many relevant scenarios.
In summary, we introduced a systematic approach to multipartite entanglement detection in many-body systems from the knowledge of the average values of arbitrary observables, whose polynomial cost at every level is guaranteed with no assumptions about the structure of the data. By benchmarking it on realistic many-body data, we are able to show that this approach has a wide range of applicability, and is able to recover and generalize several entanglement witnesses tailored to many-body systems of immediate experimental relevance.
In Section II, we present our framework for data-driven entanglement detection. In Section III we present an illustrative simple example for a Bell pair. In Section IV, we apply our method to theoretical data of realistic many-body systems, both for quench experiments, and for thermal equilibrium states. Section V displays our conclusions. More technical considerations on our method are given in Appendix A, Appendix B contains the detailed derivation of a new bipartite entanglement witness discovered through our approach, while Appendix C derives the entanglement criteria of refs. [25,26] within our framework.

II. FRAMEWORK
Some of the technical derivations of our entanglement detection method are similar to the approach presented in ref. [23]. For the sake of giving a comprehensive and self-contained description, we give complete introduction here, specialising it to the considered many-body setting. We focus on a system composed of N qubits (denoted i ∈ {1, 2, . . . , N } =: [N ]), described by an unknown quantum stateρ. We assume that the average values of several quantum observablesÔ r are known. Our ultimate goal is to prove, only from the knowledge of these average values, that the quantum stateρ is entangled.This will be achieved by exhibiting a specific entanglement witness operator, in the form of a linear combination of theÔ r operators, which is violated by the data under consideration. These average values are either obtained by directly measuring the observables in question, or are inferred from other measurements [22]. Throughout this work, entanglement is defined as the impossibility to decompose the many-body density matrix as a statistical mixture of product states over individual qubits. This encompasses the situation where all qubits are individually addressed, as is the case in typical quantum computing or quantum simulation applications; but also the situation where the qubits are two-level subspaces of indistinguishable particles, as is often the case in atomic ensemble experiments. In this latter case, where the two levels can be either two spatial modes or two internal states, correlations among the qubits are only probed via collective measurements (typically, fluctuations of collective spin observables) [4]. Such information can be naturally incorporated in our approach in order to probe entanglement among the particles.
Available quantum data. For simplicity, in what follows we assume that some one-and two-body correlations have been obtained (the more general situation, where the average value of an arbitrary collection of operators is known, is discussed in Appendix A 3). We introduce the following notations for these data: whereX i ,Ŷ i ,Ẑ i denote the qubit Pauli matrices. Notice that some of these correlators might be unknown. For instance, cross-terms such as C XY ij or C XZ ij , whose measurement require individual addressing of the qubits, are often more challenging to infer than C XX ij , C Y Y ij , C ZZ ij which can be measured via global rotations of all qubits before measuring in a fixed basis. We therefore allow for an incomplete data set Dρ = {C r } R r=1 composed of only a subset of all possible correlators [we introduce the generic notation C r := Tr(ρÔ r ) to denote either C a i or C ab ij for some 1 ≤ i < j ≤ N ; and some a, b ∈ {X, Y, Z}]. The method we develop verifies necessary conditions which are obeyed by Dρ if it can be reproduced by a separable state. The violation of any of these conditions leads our algorithm to produce a specific entanglement witness, tailored to the data under investigation, whose violation proves that the stateρ is entangled (see Fig.1 for a pictorial representation).
Sufficient conditions for entanglement. By definition, a state is separable (i.e. not entangled) if it can be decomposed as a statistical mixture of product states: Here, we have represented the local statesρ ni in a Blochsphere picture: where the local variables n i = (x i , y i , z i ) satisfy: It follows that p[{n i }] ≥ 0 can be seen as a joint probability distribution over the unit (Bloch-sphere) vectors n i for all the qubits i = 1, . . . , N . Using that: {C r } obs entanglement witness FIG. 1: Geometrical representation of the proposed entanglement detection framework. By arranging the observed data Dρ = {Cr} R r=1 as a vector, one can represent them as a point in a R-dimensional space. Among all valid quantum data (represented as the dark orange convex set), one can identify the separable set, namely the convex subset of data which can be reproduced with a separable state (yellow set). We consider an efficient way to characterize a strict superset, corresponding to the data {Cr} Γ 0 compatible with a positivesemidefinite correlation matrix Γ (light orange set). Such a set contains {Cr}sep; hence, if the data does not pass the Γ 0 test, then Dρ necessarily lies outside of the separable set, constituting a proof of entanglement. The method also provides an entanglement witness, i.e. a hyperplane separating the observed data from the separable set.
for a, b ∈ {x, y, z}. We denoted p i (resp. p ij ) the marginal distribution over the i-th qubit (resp. the (i, j) pair); and introduced the notation · · · for expectation values over the p distribution. In order to detect entanglement, one has to prove that the observed correlations {C a i , C ab ij } cannot be reproduced by the expressions (6) for any choice of joint probability distribution p[{n i }]. Crucially, one can derive conditions which are necessarily satisfied if a distribution p[{n i }] reproducing the data exists -conditions whose violation is hence sufficient to conclude that the stateρ is entangled. In order to do so, one first defines the set of classical variables m = (1, x 1 , y 1 , z 1 , . . . , x N , y N , z N ), and construct the correlation matrix Γ α,β = m α m β over the p distribution (more general choices of sets m might be considered, and are discussed in Appendix A 3). The correlation matrix Γ satisfies the following properties: • it is symmetric, Γ = Γ T .
• Some of its entries correspond to the observed data C a i and C ab ij . For instance, for m α = x i and m β = x j , then Γ α ,β = x i x j .
• Lastly, some remaining entries obey additional linear constraints, because of the condition (4). In particular, we have As an example, consider the case in which the one-body terms The corresponding Γ reads: Notice that we have marked in blue the entries replaced with the available data. All the other entries (black terms · · · ) are unknowns which represent unobserved correlations over the p distribution. If other correlators were known (e.g. C XY 12 ), they would simply replace the corresponding free variables in Eq. (7) (namely x 1 y 2 ): reducing the number of free variables makes it harder to complete the matrix Γ 0, and therefore makes it easier to detect entanglement. Notice that correlations such as x i y i have no experimental meaning in quantum physics, as they involve the simultaneous measurement of two incompatible observables, namelyX i andŶ i on the same qubit. However, they represent perfectly welldefined quantities if the state is separable, as (classical) expectation values over the p distribution. Therefore, if the state is separable, it must be possible to complete the Γ matrix with such unobserved correlations, such that Γ 0. Crucially, solving this problem is a so-called semi-definite program [24], for which efficient convexoptimization algorithms are available. As we illustrate in Section IV, the scalability of the method allows one to detect entanglement in systems of hundreds of qubits in a data-agnostic manner -that is, without a priori knowing the suitable entanglement criteria.
Construction of an entanglement witness. Importantly, if the matrix Γ 0 cannot be completed, the theory of semi-definite programming allows one to derive an entanglement witness of the form: where the sum runs only over the available data. As discussed in Appendix A 2, inequality (8) is satisfied by all separable states, while the data under investigation are such that: ultimately proving that the quantum state generating these data is entangled. The parameter λ > 0 in Eq. (8) can be interpreted as the noise robustness of the data. Indeed, if the quantum stateρ generating the data is mixed with white noise:ρ → (1 − λ)ρ + λ1/D with D = 2 N the dimension of the Hilbert space, then using the fact that Pauli observables are traceless, we have The parameter λ thus exactly quantifies the maximal amount of white noise which can be tolerated before entanglement detection becomes impossible with Eq. (7). Hence, by using an SDP to minimize the noise strength λ for which the noisy data becomes compatible with a Γ 0, one obtains the maximally robust witness possible with the method (see Appendix A 2 for details).
A converging hierarchy of conditions. The presented approach can be extended to include also higher-order correlators.
As further discussed in Appendix A 3, one may consider the set of classical where a i , b j , c k are any components of the local classical variables {n i }. One then constructs the (PSD) correlation matrix Γ α,β = m α m β over the p distribution. Verifying the compatibility of the data with Γ 0 is again a semidefinite program, which can be solved at a computational (memory) cost scaling at most as O(length(m ) 2 ). The matrix Γ [Eq. (7)] is obtained as a submatrix of Γ , and therefore the condition Γ 0 is stronger than Γ 0. Including in m all monomials up to degree l = 1, 2, 3, . . . , one defines a systematic hierarchy of positive-semidefinite conditions which are necessarily obeyed if the underlying stateρ is separable. Crucially, as further discussed in Appendix A 3, if no separable state can reproduce the data, then there exists a finite degree l such that the data fail to fulfill the corresponding condition Γ 0 -this property is a consequence of the variables {n i } being compact: |n i | = 1 for all i. Therefore, the approach presented here defines, in the limit l → ∞, a converging hierarchy of outer approximations to the set of separable states, exhausting the capability of a given data set to demonstrate multipartite entanglement. The computational cost O((3N ) 2l ) is strictly polynomial at each relaxation level. One may regard such a hierarchy as an instance of Lassere's relaxation of the moment problem for the probability distribution p[{n i }] [24,37]. Notice that in practice, the computational cost of higher-level tests (l ≥ 2) increases rapidly, especially for hundreds of qubits. However, we provide in Section IV compelling evidence of the efficiency and tightness of Eq. (7), which represents the l = 1-level of the hierarchy, to detect entanglement in many-body systems in a flexible, unbiased and scalable manner.
Invariance under partial transposition. It is interesting to notice that our criteria are independent of the partial-transposition (PT) criteria [32,33,38,39]: a stateρ is compatible with Eq. (7) if and only if (iff) the stateρ PT is compatible with Eq. (7), whereρ PT is obtained by applying partial transposition on any subset of qubits. Indeed, PT leaves invariant the Pauli matriceŝ X i andẐ i , whileŶ i is changed into −Ŷ i . Therefore, all correlations involvingŶ i are changed into their opposite. The corresponding matrix Γ PT is then obtained from Γ [Eq. (7)] by a simple change of basis, in which ( for the qubits where PT is applied. Therefore, Γ PT can be completed as a PSD matrix iff Γ can be completed as a PSD matrix. As discussed in Section A 3, this simple observation can be extended to the complete hierarchy of criteria derived via our approach.

III. SIMPLE TWO-QUBIT EXAMPLE
As a first illustration of the method, we consider an isotropic Werner state [40], namely a statistical mixture of white noise with a spin singlet: Let us show that Eq. (7) is tight for the Werner state, namely that it detects entanglement whenever λ < 2/3. The Werner state is SU (2) invariant, and one finds C a 1 = C a 2 = 0 (for a ∈ {X, Y, Z}) and C ab 12 = −δ ab (1 − λ). However, this detailed property, impossibly to exactly fulfill in an experiment, is not needed to demonstrate entanglement with our method. It turns out to be sufficient to consider only c := C XX 12 + C Y Y 12 + C ZZ 12 as available data. As discussed in Appendix A 1, if only c is known and without making any assumption about the underlying quantum state, one may symmetrize the distribution p({x i , y i , z i }) aimed at reproducing the data with a separable state, Eq. (2). This leads us to drastically simplify Eq. (7) as: Reorganizing the lines and columns, the matrix Γ is block diagonal, and is PSD iff all blocks are PSD, that is, iff 1 c c 1 0, iff |c| ≤ 1. This establishes that entanglement is detected whenever |C XX 12 + C Y Y 12 + C ZZ 12 | > 1, which is a (tight) entanglement witness already well known in the literature, and which is recovered by our approach. The witness is tight, as the product state | ↑↑ is s.t. C XX 12 + C Y Y 12 + C ZZ 12 = C ZZ 12 = 1. In the case of the Werner state, we have c = −3(1 − λ), from which we recover the known result that the state (10) is entangled for λ < 2/3. Notice that Eq. (7) is not tight for all twoqubit states: by producing random two-qubit states, we could find entangled states (as detected by the concurrence criterion [41]) which are nevertheless compatible with a PSD correlation matrix as in Eq. (7).

IV. ROBUST DETECTION OF ENTANGLEMENT IN MANY-BODY SYSTEMS
We have then chosen to benchmark our entanglement detection method on paradigmatic lattice quantum spin models. Motivated by an ultracold atoms experiment [42] realized a few years ago [43], we have first focused on the entanglement generated by a single impurity propagating along a one-dimensional XX chain (Section IV A).
As a main result, we found that the robustness of entanglement detection can be increased by about one order of magnitude using our method as compared to existing criteria, using the same data as collected in the experiment of Ref. [43]. We have then considered a twodimensional system, where entanglement is generated by a XX Hamiltonian with 1/r 3 interactions, from an initial state with all spin polarized along X (Section IV B). This example of a many-body quench with power-law interactions is especially motivated by Rydberg arrays [17,19,20,[44][45][46][47][48], ultracold magnetic atoms [49][50][51], nitrogen-vacancy centers in diamond [52,53] and trapped ions systems [10,12,[54][55][56][57][58][59], where related spin Hamiltonians have been implemented. In this case, our algorithm lead us to discover a wide family of entanglement witnesses based on components of the structure factor, which extend similar criteria reported previously in the literature, and which are especially suited to detect entanglement in out-of-equilibrium situations. Finally, we have explored the possibility to detect bipartite entanglement with our method, focusing on thermal-equilibrium states of the Heisenberg model and of the transverse-field Ising model in one dimension (Section IV C). Overall, these examples validate our method as a robust, flexible, effective and highly scalable approach to detect multipartite entanglement from partial information, as is currently collected in intermediate-scale quantum simulators and computers.
Importantly, since our main objective is to derive witnesses that can tolerate a realistic amount of noise, we always solve the entanglement SDP test as a noise robustness problem. That is, we mix the considered quantum states ρ with white noise, modelled as a completely mixed state:ρ → (1 − λ)ρ + λ1/2 N . The noise robustness λ * is then defined as the value of λ above which entanglement is not detected any more by the moment matrix criterion. By doing so, the witness obtained by the dual of the SDP tolerates, by construction, at least the amount of noise λ * (cf. Sec. II and Appendix A 2 for details). To generate the numerical SDP problems, we use the software Ncpol2sdpa [60], and we solve the SDPs with Mosek [78]. We release an open source code [79], which allows to recover results of Section IV A, and can be adapted to probe more general data.
A. Single-spin-flip in a one-dimensional chain We consider a one-dimensional ferromagnetic XX chain with nearest-neighbour interactions: with J = 1 a global energy scale, and with periodic boundary conditions. As initial state, we consider the ferromagnetic state |Ψ 0 = ⊗ i | ↑ . We assume that at time t = 0, the spin at i = 0 is flipped into | ↓ [43,61]. This central excitation then propagates along the chain under the XX Hamiltonian. AsĤ XX conserves the total magnetization along Z, the dynamics occurs in the N -dimensional manifold of states generated by the lowering operator]. Even though this simple quench is in essence a single-particle problem, multipartite entanglement is generated across the entire system. In the experiment of Ref. [43], the propagation of entanglement was observed through a lower-bound to the pairwise concurrence [61], which measures the entanglement of the two-body reduced stateρ ij [41]. Here, our main result is that using the same information as in the experiment of Ref. [43] [namely, the transverse correlations , the magnetization C Z i and the longitudinal correlations C ZZ ij ], more robust detection of entanglement is possible thanks to our method.
In order to theoretically compute the spin-spin correlations, we assume periodic boundary conditions on a chain of N = 64 spins (these choices have no visible effect on the results if the time is not long enough for the excitation to travel across the whole chain). This leads to: As discussed in Appendix A 1, in order to implement our algorithm, we may use the symmetries of the problem to drastically reduce the number of non-zero variables in Eq. (7), greatly improving the scalability. The resulting witness at each time, reconstructed via the algorithm described in Appendix A 2, is then tailored to the structure of correlations at that particular time, and follows the propagation of the excitation along the chain. The witness operator is of the formŴ Fig. 2, we plot for time tJ = 10 the correlations used as input to the SDP algorithm (upper row), and the coefficients of the corresponding entanglement witness (lower row). Both our witness and the concurrence lower-bound maximized over all pairs use the exact same data to detect entanglement. In order to compare their respective strength in a meaningful way, we have chosen to compute the noise robustness of the concurrence lower-bound as well. In Fig. 2(g), we plot the evolution as a function of time of the noise robustness for both our witness, and for the concurrence lower-bound [61] as measured in the experiment of Ref. [43]. The SDP witness is about one order of magnitude more robust againt white noise than the concurrence lower-bound. Clearly, beyond the quantitative information provided by the noise robustness, the structure of the witness also provides qualitative insight into the distribution of multipartite entanglement across the system. In particular, as is apparent in the transverse coefficients w ⊥ ij [ Fig. 2(e)], the qubits whose contribution to the witness is the largest are located close to i = −j = ±vt (with v = 1 the group velocity of the excitation). This feature is also captured by the two-body concurrence which is maximal for this pair of qubits. However, the precise contribution of other correlations is crucial to obtain a robust entanglement witness, as established by our data-agnostic approach. Finally, we notice that the separable bound as obtained from the SDP is tight, as we could always saturate this bound by a variational search over product states.
The single-particle nature of this problem is reflected in the fact that multipartite entanglement is progressively diluted throughout the system while the excitation, initially localized at i = 0, spreads across the whole chain. As a consequence, at long times, the robustness of the violation decreases to zero for large systems [ Fig. 2(g)].
In the following example, we consider instead a genuine many-body problem where the entanglement generated by the unitary dynamics is robust at all times.
B. Many-body quench dynamics in a two-dimensional power-law XX model We now consider a two-dimensional XX model with 1/r 3 interactions: where r ij denotes the distance between spins i and j, arranged over a N = L × L square lattice. We consider both open-and periodic boundary conditions with N = 400 spins. We set J = 1; and the transverse field h = 0.5 is introduced for technical reasons (see below). This model is of direct relevance both to Rydberg arrays [48], and to trapped ions [59]. As initial state, we consider a ferromagnetic state along X: For this particular initial state, the dynamics is invariant under the changê H XX → −Ĥ XX ; and |Ψ 0 represents the mean-field ground state of −Ĥ XX [62]. The dynamics is then wellapproximated by a semiclassical spin-wave approach [63], involving bosonic gaussian states, whose stability is further enhanced by introducing the symmetry-breaking term h iX i . We would like to emphasize that simulating the exact dynamics of a quantum many-body system is a central issue for all numerical approaches, and we selected this particular example, amenable to a semiclassical treatment, for the sake of illustrating the suitability of our entanglement-detection method to large-scale systems with no translation invariance, as investigated in existing experimental platforms. Ultimately, our method unveils the (in)compatibility of a given set of correlations with a separable state, and the way in which these correlations were obtained (through exact computation, using some approximations as we achieve here through a spin-wave approach, or experimentally) is totally irrelevant to the method itself. As input data, we have used the one-body terms C X i (for all qubits i) and the two-body terms C aa ij (for a = X, Y, Z and all pairs i < j). Once again, we used symmetries to reduce the number of free variables in the SDP algorithm (see Appendix A 1, and Appendix A 2 for details on the algorithm used to reconstruct an entanglement witness from the data). We have first considered systems with periodic boundary conditions, such that the data are translationally invariant (TI). In this case, we could analyze analytically the witnesses found by our algorithm, and generalize them to a whole family of entanglement witnesses in its own right. We have then considered systems with open boundary conditions, illustrating the scalability of our approach to detect entanglement in generic (non-TI) systems with hundreds of qubits in a data-agnostic manner.
A family of entanglement witnesses. Investigating TI systems, and generalizing the entanglement witnesses reconstructed by our algorithm, we found the following family of witnesses: where φ a (j) are arbitrary local phases, potentially depending on the spin direction a. The proof of the separable bound is straightforward: assuming that the state is separable, we may introduce the local variables x i , y i , z i parametrizing the local quantum states (Section II). Using that N = i (x 2 i + y 2 i + z 2 i ), we then have: The witness of Eq. (15) turns out to be very similar to existing results in the literature [64][65][66]. There are however two important differences: on the technical side, the witness of Eq. (15) involves local phases φ a (j) which may depend on the spin direction a (this possibility was not pointed out in the mentioned references [64][65][66]); and on the conceptual side, it was inferred from our algorithm in a completely data-agnostic way, as the optimal witness for TI data at the first relaxation level of our hierarchy. In this case, the local phases are of the form φ a (j) = k a · r j with r j the position of the j-th subsystem. This leads to the structure factor S a k = N −1 j,j e ik·(r j −rj ) C aa jj (notice that we have included the j = j term in the summation, corresponding to a term C aa jj = 1). Although we discovered these witnesses focusing on two-dimensional systems, they can be naturally extended to arbitrary geometries. In terms of components of the structure factors, the entanglement witness of Eq. (15) reads: Notice that we have defined the structure factors in terms of qubit observables, which are twice the spin observables typically used in condensed-matter physics; this leads to a factor 4 in the definition of the structure factors. Crucially, the wavevector k may be different for the X, Y and Z components of the spins. Obviously, in order to detect entanglement, it is optimal to choose, for each spin component a, the direction k a where the structure factor is minimal, leading to: It is in this form that the witnesses have been discovered via our algorithm, investigating the correlations generated by the dynamics in the power-law XX model. We have then generalized this result to Eq. (15). The witness Eq. (17) is then violated if the fluctuations at these optimal wavevectors are suppressed below the separable bound 2. Physically, these witnesses detect a generalized form of spin squeezing, especially suited to the many-body systems where different components of the structure factor can be measured. In several quantum simulators, the structure factors are reconstructed by Fourier transform of the real-space correlations. In condensed-matter systems, correlations are typically measured directly in momentum space via neutron scattering. During the out-of-equilibrium dynamics, these optimal wavevectors may vary over time in different ways for different spin components; and therefore the witnesses of Eq. (15)-(16)-(17) offer a large flexibility to detect entanglement, independently of the specific SDP algorithm we used to discover them. Finally, we notice that the witness of Eq. (17) can be extended to spin-s systems. We define in general the structure factor as S a k = N −1 j,j e ik·(r j −rj ) Tr[ρŜ a jŜ a j ] withŜ a j the spin-s observable in direction a for subsystem j. Using that This generalizes a result of Ref. [67] to arbitrary components of the structure factors.
Numerical results. In Fig. 3(a), we plot the structure factor at time tJ = 1, for wavevectors k = (k, k), in a 20 × 20 square lattice with open boundary conditions. The minimal value of the structure factor for each spin component (marked by a circle on the figure) then enters the TI witness of Eq. (17). In Fig. 3(b), we plot the noise robustness of the (non-TI) witness found by our algorithm as a function of time.  17), whose reconstruction is illustrated on panel (a) for time t = 1/J. evaluated at the (time-dependent) optimal wave-vectors (k X , k Y , k Z ). While the TI witness reaches a noise robustness of about 0.5, the (non-TI) optimal witness tailored to the (non-TI) correlations reaches more than 0.8 noise robustness. However, we could not find an analytical expression for these non-TI data-driven witnesses. By a variational search over separable states, we could however verify that the separable bound obtained by our SDP algorithm was always tight.

C. Bipartite witnesses
Finally, we show that the very same SDP algorithm outlined in Sec II can be adapted with no additional computation cost to detect bipartite entanglement along any splitting of the system in two parts. We therefore consider a partition of the N qubits into two halves A and B, and as input data, we consider single-qubit terms C a i and only inter-AB correlations C ab ij where i ∈ A and j ∈ B. It is straightforward to notice that if non-fullseparability can be proved from this knowledge, then the state must be bipartite entangled. Indeed, if the AB state is bipartite separable:ρ bisep ]. One can verify that ifρ bisep AB reproduces the data, so doesρ fullsep AB ; conversely, proving non-full-separability from these data implies bipartite entanglement.
Heisenberg model and transverse-field Ising models. We use the above idea to investigate bipartite entanglement in the transverse-field Ising and Heisenberg chain at finite temperature, for N = 64 spins. The Heisenbeg chain is described by the Hamiltonian: (18) and transverse-field Ising chain by: In Eqs. (18) and (19), we have assumed periodic boundary conditions. The parameter J is an overall energy scale (set to J = 1 in our computations), and g in Eq. A family of bipartite entanglement witnesses.
Having first considered a partition of the form AAA . . . AA|BBB . . . BB (namely, A = {0, 1, 2, . . . , N/2 − 1} and B = {N/2, . . . , N − 1}), we have noticed that entanglement was detected iff the nearest-neighbour two-body reduced density matrix ρ N/2−1,N/2 was itself entangled (as detected by the concurrence [41]). While illustrating the relatively shortrange nature of entanglement in these thermal states, we could not go beyond the mere witnessing of entanglement among nearest-neighbours. We have therefore considered a partitionning maximizing the AB interface, that is: where φ a (j) are arbitrary local phases, and with coeffi- cients given by: Notice that for i ∈ A and j ∈ B, r = j − i is an odd integer, in which case we have the simplified expression K r = 2(−1) (r−1)/2 N tan(πr/N ) . The violated witnesses are then of the form The proof of the witness inequality (24) is given in Appendix B. For the Heisenberg chain, the optimal choice of the phases is φ a (i) = 0. For the Ising chain, it is φ a (i) = 0 for i ∈ A, and φ X (j) = φ Z (j) = π for j ∈ B, and φ Y (j) = 0. As illustrated in Fig. 4, the bipartite entanglement witnesses of Eq. (24) allow one to detect entanglement in regimes where all two-body reduced density matrices are separable (as measured by the concurrence [41]).
Comparison with the structure factor witnesses. For the sake of completeness, we have also evaluated the multipartite entanglement witness based on the structure factor [Eq. (17)] discussed in Section IV B. As this witness involves also intra-A and intra-B correlations, and is based on the same SDP criterion, it detects entanglement at strictly higher temperatures than the bipartite witness of Eq. (24).
The optimal wavevectors k X , k Y and k Z are found from the following observations. The Heisenberg model develops antiferromagnetic correlations at low temperature, leading to a peak at k = π in the structure factor. Concomitantly, fluctuations of the uniform magnetization (at k = 0) are suppressed, reaching a spinsinglet state in the ground state (S a 0 = 0). At all temperatures, the structure factor is always minimal at k X = k Y = k Z = 0, and in this case the optimal witness is simply a S a 0 ≥ 2. As previously reported in Ref. [35], it is violated up to a temperature T /J ≈ 1.4 [ Fig. 4(a)]. In contrast, at low temperature and around the quantum critical point g = 1 [68], the Ising model develops ferromagnetic correlations for the Z component of the spin, and the structure factor is minimal at k Y = π. The uniform magnetization along Y is slightly squeezed below the standard quantum limit [69], and the structure factor is minimal at k Y = 0. Finally, correlations in the X direction (the direction of the transversefield) are strongly ferromagnetic, and are suppressed at k X = π. We have found that these choices are optimal throughout the phase diagram. As illustrated in Fig. 4(b), above the quantum critical point g = 1 the optimal witness S X π + S Y 0 + S Z π ≥ 2 is violated for temperatures T /J 0.38. In contrast a criterion based on the quantum Fisher information (based on the dynamical structure factor for ZZ correlations [70], and which is considerably more challenging to estimate, both in theory and in experiments [71,72]) is violated for T /J 0.11; and a witness based on the same data as used in the present work [35], and optimized for T /J = 0.28, is violated for T /J 0.31. The entanglememt witness of Eq. (17) discovered in the present work offers therefore both a simpler and more robust criterion for the detection of multipartite entanglement in many-body systems, as compared to existing criteria proposed so far in the literature.

V. CONCLUSION
We have implemented a systematic, scalable and flexible approach to detect multipartite entanglement in many-body systems. Assuming that the knowledge of some average values of many-body observables are known, one can build a correlation matrix whose entries reproduce these data. Under the assumption that the state is separable, positive semidefinite constraints must be obeyed by the correlation matrix. Verifying these constraints is efficiently achieved via semidefinite programming techniques, yielding a data-tailored entanglement witness violated by the data under consideration. We have illustrated the scalability of this approach in some paradigmatic examples of many-body systems, demonstrating for instance how our approach can easily deal with systems of hundreds of qubits when two-body correlations are used as input data.
By choosing to perform the proposed entanglement test as a noise robustness problem, we were able to show that the corresponding entanglement witnesses can tolerate realistic amount of noise in many physically relevant many-body scenarios. We have also shown that the specific entanglement witnesses discovered via our approach can sometimes be analyzed analytically, leading to explicit entanglement witness of independent relevance [see e.g. Eq. (15), (17) or (24)].
Within our framework, one can probe the (in)compatibility of a given set of average values with a separable state for any fixed partitioning of the system. Here, we focused mostly on a partitioning into N individual qubits, or a bipartition into two halves of N/2 qubits. Our approach can be used to detect entanglement among individually addressed qubits, as extensively demonstrated throughout the paper. It can also be used to detect entanglement among indistinguishable particles, probed e.g. via collective spin observables; as a matter of fact, all generalized spin-squeezing inequalities [26] typically used to detect entanglement in this framework are recovered as a special case by our method. It is a priori unclear if our approach can be extended to encompass also statistical mixtures of different partitionings, as required to quantify the so-called entanglement depth, or width [73].
Several research directions are now open to future works. Fist, the considered approach can be extended to qudit systems in a straightforward manner, which provides a scalable technique for the investigation of multipartite entanglement in large ensembles of qudits. Then, although we have illustrated the method by assuming that one-and two-body correlations are known, one could naturally include the knowledge of any k-point function. Including such higher-order correlations would certainly enhance the capability to detect entanglement in the examples we have presented. Studying topological phases [74], where entanglement could be revealed by stringorder-parameters, via specific entanglement witnesses inferred by our data-driven method, represents also an exciting avenue for future works. Finally, implementing our algorithm using experimental data as input would probably reveal unforeseen features of many-body entanglement. This last possibility is especially relevant in the context of quantum computation and quantum simulation, operating beyond the capabilities of classical computers [14][15][16][17][18][19][20], and where entanglement is commonly assumed to be an essential resource [2].
after the completion of the first version of this paper. This work is supported by the ERC AdG CERQUTE, the AXA Chair in Quantum Information Science, the Government of Spain (FIS2020-TRANQI, Severo Ochoa CEX2019-000910-S and Retos QuSpin), Fundacio Cellex, Fundacio Mir-Puig, Generalitat de Catalunya (CERCA, AGAUR SGR 1381 and QuantumCAT), the Austrian Science Fund (FWF) through Project number 414325145 within SFB F7104 and the Alexander von Humboldt Foundation, the Agence Nationale de la Recherche (ANR) Research Collaborative Project Qu-DICE (ANR-PRC-CES47), the John Templeton Foundation (Grant No. 61835). the data. Importantly, p 4 is such that x i p4 = 0 (since under p 4 , x i has the same probability as −x i ). Similarly, y i p4 = 0, and two-body terms are such that x i y j p4 = x i y j p4 = y i z j p4 = 0. Therefore, without loss of generality, one may impose in the SDP algorithm that x i = y i = x i y j = x i z j = y i z j = 0 for all i, j. One thus reduces the number of free variables in Eq. (7) from order O(N 2 ) to 2N (for instance, all x 2 i and y 2 i for i ∈ [N ]).
2) Let us then consider a situation where the data consist of as was the case in the study of the quench in the one-dimensional XX chain. Let us assume that a probability distribution p({x i , y i , z i }) reproduces the data.
We consider then the distribution p = p({y i , x i , z i }). As x i x j + y i y j p = x i x j + y i y j p , this distribution also reproduces the data. Hence, we may as well consider q = (p + p )/2, which also reproduces the data. As q is such that x i x j q = y i y j q for all pairs (i, j), we may assume C XX ij = C Y Y ij in the SDP algorithm without loss of generality, and without assuming that this symmetry is actually present in the experiment. Following similar arguments, the distribution p = p({−x i , y i , z i }) reproduces the data, since x i x j p = x i x j p . Considering (p + p )/2, we may impose x i = 0, and also x i y j = x i z j = 0 for all pairs (i, j). The same reasoning apply interchanging the role of x and y, leading to y i = y i z j = 0. In conclusion, we may solve the SDP in the form of Eq. (7) imposing that all terms are zero, except the data 3) Finally, if one only knows c ij := C XX ij +C Y Y ij +C ZZ ij , by a similar reasoning one may consider a distribution p({x i , y i , z i }) which is invariant under all rotations, satisfying x i = y i = z i = 0, x i y j = x i z j = y i z j = 0, and x i x j = y i y i = z i z j . In this case, there is no free variable at all in Eq. (7) [all two-body diagonal terms are 1/3, and the only non-zero off-diagonal entries ij , obtained as c ij /3]. As was illustrated in the case of the Werner state, testing entanglement via Eq. (7) simply consists of checking positive semidefiniteness of the matrix M , defined by M ij = c ij if i = j, and M ii = 1. We emphasize that this does not assume any spatial symmetry, and applies in particular to data with no translation invariance.

Robust entanglement witness from one-and two-body correlations on N qubits
Here we provide details on the way to extract an entanglement witness of the form (8) via semidefinite programming. Specifically, we show how to find the noise robustness of the entanglement contained in the data, and derive an entanglement witness via semidefinite programming. We start considering a situation where N qubits are measured, and some of the one-body terms C a i and two-body terms C ab ij have been measured. We denote generically these data as {C α } R α=1 := {C a i , C ab ij }, where α labels both the sites and measurements. We assume that in total the dataset contains R single-and two-body correlations. We define the noise robustness as the minimal value of λ such that {(1 − λ)C α } is compatible with a positive semidefinite (PSD) correlation matrix, as explained in the main text [Section II and in particular Eq. (7)]. That is, we aim at solving the problem: The (symmetric) matrix Γ = (Γ n,m ) 0≤n,m≤3N is the correlation matrix for the variables parametrizing separable states of N qubits (Section II), and is therefore PSD [Eq. (A1a)]. In Eq. (A1b), we constrain the relevant entries of the Γ matrix to reproduce the data (with a 1 − λ noise prefactor); we have introduced (n, m)(α) to denote the pair of indices (n, m) containing the data C α . Specifically [see Eq. (7)], C a i is contained in Γ 0,3(i−1)+a , and C ab ij in Γ 3(i−1)+a,3(j−1)+b . Finally, Eq. (A1c) enforces x 2 i + y 2 i + z 2 i = 1 in the Γ matrix; as this condition descends from properties of the Pauli matrices (Section II), we call it the Pauli constraint.
Standard primal form. We now rewrite in its socalled standard primal form [24] the semidefinite program (SDP) defined in Eq. (A1): We have introduced the matrix scalar product X, Y = Tr(X T Y ) = ij X ij Y ij . The involved matrices have the following block-diagonal form: where the Σ matrices are given by: The SDP (A2) is equivalent to the problem (A1) of finding the minimal noise λ for which the Γ matrix in (7) becomes positive semidefinite, for a set of noisy data Therefore, an optimal solution λ * > 0 implies that the given data are not compatible with a separable state, hence resulting in entanglement detection.
Dual form. If that is the case, one can derive an entanglement witness by considering the so-called dual problem [24] corresponding to Eq. (A2): Using the expressions of the matrices M and A's [Eqs. (A3) and (A4)], we rewrite Eq. (A5) as: We denote { w * , ( w Pauli ) * } the optimal solution to this dual problem.
Strong duality. As for all semidefinite progams [24], any primal feasible X [that is, a PSD matrix X 0 satisfying the constraints in Eq. (A2)] yields an upper bound to any dual feasible [that is, any {w α , w Pauli i } satisfying the constraints in Eq. (A5)]. Indeed, for any such feasible X and w, we have: ), X , where we used the constraints in Eq. (A2).
Using then the constraint in Eq. (A5), we notice that As X 0, and using the fact that Y, X ≥ 0 for any two PSD matrices Y and X, we conclude that M, X − α w α C α − i w Pauli i ≥ 0. In particular, the primal optimum upper bounds the dual optimum: λ * ≥ w * · C + i (w Pauli i ) * , a property known as weak duality [24]. In our case, the primal and dual optima are actually equal, as a consequence of strong duality which holds for our problem. A sufficient condition for strong duality to hold [24] is that both primal and dual problems are strictly feasible. The primal problem is stricly feasible if one can find a positive definite matrix X 0 satisfying the constraints in Eq. (A2): such strictly feasible X is readily obtained by choosing λ = 1 and Γ = 1/3 in Eq. (A3) [or, equivalently, in Eq. (A1)]. The dual problem is stricly feasible if one can exhibit some {w α , w Pauli i } satisfying the constraints in Eq. (A6) as strict inequalities: such stricly feasible w's are readily obtained as w α = 0 and w Pauli . Two important properties follow from strong duality: 1) ), X * = 0. Using the expression of M and the A's matrices [Eq. (A3)], this implies in particular (1 − w * · C)λ * = 0. Whenever entanglement is detected (λ * > 0), we have therefore w * · C = 1, so that − i (w Pauli Entanglement witness. The coefficients w * α define an entanglement witness whose separable bound is given by which holds for all separable data {C sep α } (namely, data compatible with a separable state), and which is violated by the data under consideration ( α w * α C α = 1). In order to prove this fact, we first observe that w * · C sep ≤ 1 for all separable data C sep (indeed: w * · 0 = 0 for the separable data 0, while w * · C = 1 for the non-separable data C; we conclude using the convexity of the set of separable data). Therefore, { w * , (w Pauli i ) * } represent a dual feasible [Eq. (A6)] for any separable data C sep . As a consequence of duality, this provides a lower bound to the primal optimum, which is (λ * ) sep = 0 for separable data. Hence, systematic hierarchy of such choices, converging towards the exact separable set. The l-th relaxation level is defined as: Throughout this work, we considered only the first relaxation level defined by Eq. (A9). Crucially, the corresponding hierarchy converges, in the limit l → ∞, towards the separable set. A way to see that is to interpret the hierarchy as Lassere's series of relaxation for the moment problem associated to the variables {(x i , y i , z i )}.
Since the variables satisfy the quadratic constraint (4), the relaxation meets the Archimedean condition, which is enough to guarantee convergence of the hierarchy [37]. It follows that the set of correlations that can be recovered as moments of an overall distribution p({n i }), i.e. the set of separable correlations, is obtained as the asymptotic limit of the hierarchy defined above. Therefore, if the data are incompatible with a separable state, they will be detected as entangled at a finite level of the hierarchy -althought the computation cost of high level tests quickly increases with l, as one needs to manipulate a correlation matrix of size ∼ (4N ) l in the semidefinite program.
Lastly, notice that one can straightforwardly define some hybrid levels of the hierarchy, where the entries of the vector v (l) are complemented by monomials of order higher than l. Such hybrid conditions have the flexibility of including the knowledge of a finite amount of higher-order correlations, while retaining the scalability of the computational cost given by the fixed level l.
Invariance of the hierarchy under partial transposition. In Section II of the main text, we already noticed that the level-1 relaxation is left invariant by the partial transposition (PT) of any subsystem. The key observation was that partial transposition simply amounts to a change of basis for the Γ matrix -hence, positivity of Γ is left unchanged under PT. The same observation carries over to arbitrary relaxation levels [Eq. (A12)]. Indeed, the Pauli matrices are either symmetric or antisymmetric; therefore, under PT they are either left invariant, or transformed into their opposite. In terms of the correlations of the n i variables [cf. Eq. (7)], or in terms of the choice of monomials at a given relaxation level [cf. Eq. (A12)], it simply amounts to change the corresponding variables into their opposite, which is again achieved by a change of basis. Our approach therefore represents a hierachy of conditions which are completely independent of the PPT-based criteria [32,33,38,39]. Here we prove the validity of (24) as a witness of bipartite entanglement according to an even-odd partitioning of the system. For the sake of completeness, we first revise the setting and the witness expression. We consider the bipartition A = {0, 2, 4, . . . , N − 2} and B = {1, 3, 5, . . . , N − 1}. We introduce local phases φ a (i) for a ∈ {X, Y, Z} and i ∈ [N ], and define: The coefficients K r are given by: The violated witnesses are then of the form: Proof. The proof of this inequality is as follows. Assuming that the state is fully separable, we have: We then observe that for j, j ∈ A, j − j is an even integer, and therefore, from Eq. (B3), K j−j = 1 − 2/N if j = j and K j−j = −(2/N )(−1) (j−j )/2 if j = j . Therefore, we have: By the same argument, observing the j − j is an even integer for j, j ∈ B, we have: j,j ∈B K j−j a j a j e i[φa(j)−φa(j )] ≤ r∈B a 2 r . (B8) We therefore have: Combining these inequalities for a ∈ {X, Y, Z}, we conclude that: where in the last step we have used that the constraint (4) is obeyed by fully separable states. This achieves the proof that Eq. (B4) is an entanglement witness. Combining with the observation made at the beginning of Sec. IV C, we conclude that a violation of such a witness directly implies bipartite entanglement, since its expression involves only cross-correlations between A and B subsystems.
Appendix C: Relation to previous entanglement criteria Here we show how one can recover some previously known entanglement criteria as a consequence of the PSD condition introduced in Sec. II.

Recovering the covariance matrix criterion
In this subsection, we show how the so-called covariance matrix criterion (CMC) [25] for detecting entanglement in a multiqubit state is a consequence of our approach. In order to state the CMC within the notations of this paper, we consider three qubits (the N -qubit case requires a trivial generalization). Applying the CMC requires the knowledge of all one-body and two-body correlations for all pairs, and all Pauli matrices. We denote C i the vector (C X i , C Y i , C Z i ), and C ij the 3 × 3 matrix with entries C ab ij [see Eq.
(1) for the definition of C a i and C ab ij ]. The CMC states that if the three-qubit state is fully separable, then there exist three real symmetric 3×3 matrices ρ i 0, with Tr(ρ i ) = 1, such that: On the other hand, the first level of our hierarchy [see Eq. (7)] implies that: denotes an average over the (classical) p distribution defining a separable state (Section II). Clearly, σ i 0 is a symmetric matrix of unit trace. Furthermore, positivity of Γ implies Eq. (C1). Notice that our approach is more general than the approach underlying the CMC, for at least three reasons: 1) we can naturally deal with missing entries in the correlation matrix (and actually leverage on it to strongly reduce the computational cost, see Appendix A 1); 2) the criterion of Eq. (7) which implies the CMC is only the first level of a systematic hierarchy converging towards the separable set; 3) we provide a systematic approach to incorporate the knowledge of any correlation function, beyond one-and two-body considered in the CMC.

Recovering the generalized spin-squeezing inequalities
Here, we show how the generalized spin-squeezing inequalities derived in Ref. [26] can be recovered within our approach. These inequalities are entanglement witnesses invariant under all permutations of the qubits, and may be defined in terms of first-and secondmoments averaged over all permuatations: {m a , C aa ; a ∈ {x, y, z}} with m a = N −1 N i=1 C a i and C aa = [N (N − 1)] −1 i =j C aa ij . They consist of the following eight inequalities, valid for all fully-separable states [Eq. (50) of Ref. [26]]: x + m 2 y + m 2 z ) − (N − 1)(C xx + C yy + C zz ) ≤ 1 In order to prove these inequalities within our approach, we assume that there exists a probability distribution p[{x i , y i , z i }] reproducing the data. Following the symmetrization procedure described in Appendix A 1, without loss of generality we may choose the p distribution invariant under all permutations. We then consider as monomials u = (1, x 1 , N i=1 x i ), and build the corresponding correlation matrix Γ x = u T u , where the average is over the p distribution. Using the invariance of p under permutations of the qubits, we obtain: Positivity of Γ x implies that: From these inequalities, we obtain the conditions: For the same reason, we also have: Combining these inequalities, and using the property x 2 1 + y 2 1 + z 2 1 ≤ 1, we recover the eight inequalities of Ref. [26]. The approach presented in this paper is however much more general. Notice also that a straightforward extension to qudits allows one to recover the results of Ref. [67].