# Unveiling Quantum Entanglement in Many-Body
Systems from Partial Information

Irénée Frérot, Flavio Baccari, Antonio Acín

► To cite this version:

Irénée Frérot, Flavio Baccari, Antonio Acín. Unveiling Quantum Entanglement in Many-Body Systems from Partial Information. PRX Quantum, 2022, 3 (1), pp.010342. 10.1103/PRXQuantum.3.010342. hal-03610050

HAL Id: hal-03610050

https://hal.science/hal-03610050v1

Submitted on 15 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Unveiling Quantum Entanglement in Many-Body Systems from Partial Information

Irénée Frérot,¹,²,³,* Flavio Baccari,²,† and Antonio Acín¹,⁴

¹ICFO - Institut de Ciencies Fotoniques, The Barcelona Institute of Science and Technology, Castelldefels (Barcelona) 08860, Spain

²Max-Planck-Institut für Quantenoptik, Hans-Kopfermann-Straße 1, Garching 85748, Germany

³Université Grenoble Alpes, CNRS, Grenoble INP, Institut Néel, Grenoble 38000, France

⁴ICREA - Institucio Catalana de Recerca i Estudis Avançats, Lluis Companys 23, Barcelona 08010, Spain

<figure><img src="image_2.png" /></figure>

(Received 11 August 2021; revised 2 December 2021; accepted 22 February 2022; published 15 March 2022)

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³ interactions; and thermal equilibrium states of Heisenberg and transverse-field Ising chains. In all cases, these investigations have led 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.

DOI: 10.1103/PRXQuantum.3.010342

## 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 many-body 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 [5], the ability to manipulate quantum many-body superpositions arguably represents a major endeavor 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 [6–14]. On the other hand, a growing number of experiments operate in regimes inaccessible to the best available classical simulations [15–21]—another hallmark pointing towards a genuine quantum advantage. It is commonly assumed that the intractability of classical simulations originates in the large-scale 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

Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.

<footnote>*irenee.frerot@neel.cnrs.fr</footnote>

<footnote>†flavio.baccari@mpq.mpg.de</footnote>
witnesses, for the structure of entanglement in the sys-
tem, 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 infor-
mation would require a number of measurements growing
exponentially with the system size [22,23].

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 val-
ues of (a scalable number of) arbitrary observables are
known. Starting from the same insight as in Ref. [24],
namely the connection between the compatibility of data
with (quantum) separable states and the (classical) trun-
cated 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 Ref. [24] is that if the under-
lying quantum state is separable, then the available data are
obtained as entries of a correlation matrix, which satisfies
certain positive-semidefinite constraints. Such compati-
bility conditions can be efficiently verified via so-called
semidefinite-programming (SDP) techniques [25], allow-
ing the study of systems of hundreds of qubits. The failure
for the data to pass this test serves directly as an entan-
glement 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 rel-
evant, and, in fact, is not discussed in this work. As an
illustration, we consider one- and two-body correlations
for *N* qubits, but our scheme is flexible and can incorpo-
rate 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 sys-
tem. As a matter of fact, for several of the examples we
consider, we could analytically characterize the witnesses
obtained numerically. This led us to extend some known
entanglement criteria of the literature and derive com-
pletely new ones as well. Analogously to the hierarchy
introduced in Ref. [24], the scheme we propose can be gen-
eralized 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 hierarchy—in this sense, the hierarchy
is complete.

**Comparison to previous works.** A large body of liter-
ature has already considered the problem of entanglement

detection from partial information. Some of these results
are recovered as special cases of the approach imple-
mented in this paper; some others lack the scalability
required to apply them to many-body systems; and some
alternative scalable schemes either lack the flexibility of
the present approach, or can be inconclusive. In partic-
ular, the so-called covariance matrix criterion [26] and
the generalized spin-squeezing inequalities [27], 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 [28–30]. 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 [26–30] provide
only sufficient conditions for entanglement, and there-
fore 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 incor-
porate partial knowledge about the quantum state, was
proposed based on the solution of so-called separabil-
ity eigenvalue equations [31,32]; 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 devel-
oped. Such approaches allow one to test bipartite entan-
glement criteria based on Rényi entropies [13,23], and
partial-transposition- (PT) based criteria [33,34]. How-
ever, in addition to the very high experimental require-
ments 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 infor-
mation have been developed. An entanglement-detection
method from the knowledge of two-body reduced den-
sity matrices was developed in Ref. [35], with a simi-
lar 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 aver-
age value of many-body observables. The approach of
the present paper is complementary to Ref. [36], where
the problem is solved through a mapping onto an inverse
problem of classical statistical physics, offering a sys-
tematic and scalable solution; however, this approach
could be inconclusive for particular data sets. In con-
trast, here we solve a relaxation to this problem with
an efficiency that 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 Ref. [37]. However, in contrast to ours,
the method in Ref. [37] 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. [24]. What we develop here can be seen as a complementary separability test, based on the same conceptual premises, but with a different physical motivation. While Ref. [24] aims at constructing—if it exists—a 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 Ref. [24] 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 Ref. [24], it has no guarantee to be robust against noise. Lastly, we notice that the element of randomness in the objective function in Ref. [24] 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 Sec. II, we present our framework for data-driven entanglement detection. In Sec. III we present an illustrative simple example for a Bell pair. In Sec. 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. [26,27] within our framework.

## II. FRAMEWORK

Some of the technical derivations of our entanglement detection method are similar to the approach presented in Ref. [24]. For the sake of giving a comprehensive and self-contained description, we give a complete introduction here, specializing it to the considered many-body setting. We focus on a system composed of $N$ qubits (denoted $i \in \{1, 2, \dots, N\} =: [N]$), described by an unknown quantum state $\hat{\rho}$. We assume that the average values of several quantum observables $\hat{O}_r$ are known. Our ultimate goal is to prove, only from the knowledge of these average values, that the quantum state $\hat{\rho}$ is entangled. This will be achieved by exhibiting a specific entanglement witness operator, in the form of a linear combination of the $\hat{O}_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 [23]. 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:

$$ \begin{cases} C_i^X = \mathrm{Tr}[\hat{\rho}\hat{X}_i] \\ C_i^Y = \mathrm{Tr}[\hat{\rho}\hat{Y}_i] \\ C_i^Z = \mathrm{Tr}[\hat{\rho}\hat{Z}_i] \end{cases} \quad \begin{cases} C_{ij}^{XX} = \mathrm{Tr}[\hat{\rho}\hat{X}_i\hat{X}_j] \\ C_{ij}^{XY} = \mathrm{Tr}[\hat{\rho}\hat{X}_i\hat{Y}_j] \\ \vdots \end{cases}, \quad (1) $$

where $\hat{X}_i$, $\hat{Y}_i$, $\hat{Z}_i$ denote the qubit Pauli matrices. Notice that some of these correlators might be unknown. For instance, cross terms such as $C_{ij}^{XY}$ or $C_{ij}^{XZ}$, whose measurement require individual addressing of the qubits, are often more challenging to infer than $C_{ij}^{XX}$, $C_{ij}^{YY}$, $C_{ij}^{ZZ}$, which can be measured via global rotations of all qubits before measuring in a fixed basis. We therefore allow for an incomplete data set $\mathcal{D}_{\hat{\rho}} = \{C_r\}_{r=1}^R$ composed of only a subset of all possible correlators [we introduce the generic notation $C_r := \mathrm{Tr}(\hat{\rho}\hat{O}_r)$ to denote either $C_i^a$ or $C_{ij}^{ab}$ for some $1 \le i < j \le$
N; and some $a, b \in \{X, Y, Z\}$. The method we develop verifies necessary conditions, which are obeyed by $D_{\hat{\rho}}$ 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 $\hat{\rho}$ 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:

$$ \hat{\rho}_{\text{sep}} = \prod_{i=1}^{N} \int_{|\mathbf{n}_i|=1} p[\{\mathbf{n}_i\}] \otimes_{i=1}^{N} \hat{\rho}_{\mathbf{n}_i}. \quad (2) $$

Here, we represent the local states $\hat{\rho}_{\mathbf{n}_i}$ in a Bloch-sphere picture:

$$ \hat{\rho}_{\mathbf{n}_i} = |\mathbf{n}_i\rangle\langle\mathbf{n}_i| = \frac{1}{2}(1 + x_i\hat{X}_i + y_i\hat{Y}_i + z_i\hat{Z}_i), \quad (3) $$

where the local variables $\mathbf{n}_i = (x_i, y_i, z_i)$ satisfy

$$ \forall i \in [N] \quad \mathbf{n}_i^2 = x_i^2 + y_i^2 + z_i^2 = 1. \quad (4) $$

It follows that $p[\{\mathbf{n}_i\}] \ge 0$ can be seen as a joint probability

<figure><img src="image_3.png" /><figcaption>FIG. 1. Geometrical representation of the proposed entanglement detection framework. By arranging the observed data $D_{\hat{\rho}} = \{C_r\}_{r=1}^R$ 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 $\{C_r\}_{r \ge 0}$ compatible with a positive-semidefinite correlation matrix $\Gamma$ (light orange set). Such a set contains $\{C_r\}_{\text{sep}}$; hence, if the data does not pass the $\Gamma \ge 0$ test, then $D_{\hat{\rho}}$ 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.</figcaption></figure>

distribution over the unit (Bloch-sphere) vectors $\mathbf{n}_i$ for all the qubits $i = 1, \dots, N$. Using that

$$ \left\{ \begin{aligned} x_i &= \operatorname{Tr}[\hat{\rho}_{\mathbf{n}_i} \hat{X}_i] \\ y_i &= \operatorname{Tr}[\hat{\rho}_{\mathbf{n}_i} \hat{Y}_i] \\ z_i &= \operatorname{Tr}[\hat{\rho}_{\mathbf{n}_i} \hat{Z}_i] \end{aligned} \right. , \quad (5) $$

one may express the correlators [Eq. (1)] in a separable state [Eq. (2)] as classical expectation values over the p distribution:

$$ C_i^A = \int_{|\mathbf{n}_i|=1} p_i(\mathbf{n}_i) a_i := \langle a_i \rangle, \quad (6a) $$

$$ C_{ij}^{AB} = \int_{|\mathbf{n}_i|=1} \int_{|\mathbf{n}_j|=1} p_{ij}(\mathbf{n}_i, \mathbf{n}_j) a_i b_j := \langle a_i b_j \rangle \quad (6b) $$

for $a, b \in \{x, y, z\}$. We denoted $p_i$ (respectively, $p_{ij}$) the marginal distribution over the $i$th qubit [respectively, the $(i, j)$ pair]; and introduced the notation $\langle \cdots \rangle$ for expectation values over the p distribution.

In order to detect entanglement, one has to prove that the observed correlations $\{C_i^a, C_{ij}^{ab}\}$ cannot be reproduced by the expressions (6) for any choice of joint probability distribution $p[\{\mathbf{n}_i\}]$. Crucially, one can derive conditions, which are necessarily satisfied if a distribution $p[\{\mathbf{n}_i\}]$ reproducing the data exists—conditions whose violation is hence sufficient to conclude that the state $\hat{\rho}$ is entangled. In order to do so, one first defines the set of classical variables $\mathbf{m} = (x_1, y_1, z_1, \dots, x_N, y_N, z_N)$, and construct the correlation matrix $\Gamma_{\alpha, \beta} = \langle m_{\alpha} m_{\beta} \rangle$ over the p distribution (more general choices of sets $\mathbf{m}$ might be considered, and are discussed in Appendix A 3).

The correlation matrix $\Gamma$ satisfies the following properties:

(a) It is symmetric, $\Gamma = \Gamma^T$.

(b) It is positive semidefinite (PSD, i.e., $\Gamma \ge 0$) by construction: indeed, for any vector $\mathbf{v}$, we have $\mathbf{v}^T \Gamma \mathbf{v} = \langle (\sum_{\alpha} m_{\alpha} v_{\alpha})^2 \rangle \ge 0$.

(c) Some of its entries correspond to the observed data $C_i^a$ and $C_{ij}^{ab}$. For instance, for $m_{\alpha'} = x_i$ and $m_{\beta'} = x_j$, then $\Gamma_{\alpha', \beta'} = \langle x_i x_j \rangle$.

(d) Lastly, some remaining entries obey additional linear constraints, because of the condition (4). In particular, we have $\langle z_i^2 \rangle = 1 - \langle x_i^2 \rangle - \langle y_i^2 \rangle$ for all i.

As an example, consider the case in which the one-body terms $C_i^X, C_i^Y, C_i^Z$ have been measured, together with two-body terms $C_{ij}^{XX}, C_{ij}^{YY}, C_{ij}^{ZZ}$. The corresponding $\Gamma$ reads
$$ \Gamma = \Gamma^T = \begin{pmatrix} 1 & C_1^X & C_1^Y & C_1^Z & C_2^X & C_2^Y & \dots & C_N^X & C_N^Y & C_N^Z \\ \langle x_1^2 \rangle & \langle x_1 y_1 \rangle & \langle x_1 z_1 \rangle & C_{12}^{XX} & \langle x_1 y_2 \rangle & \dots & C_{1N}^{XX} & \langle x_1 y_N \rangle & \langle x_1 z_N \rangle \\ \langle y_1^2 \rangle & \langle y_1 y_2 \rangle & \langle y_1 z_2 \rangle & C_{12}^{YY} & \langle y_1 x_2 \rangle & \dots & \langle y_1 x_N \rangle & C_{1N}^{YY} & \langle y_1 z_N \rangle \\ \dots & \dots & \dots & \langle x_1^2 \rangle & \langle x_2 y_2 \rangle & \dots & C_{2N}^{XX} & \langle x_2 y_N \rangle & C_{2N}^{ZZ} \\ \dots & \dots & \dots & \langle y_1^2 \rangle & \langle y_2 y_2 \rangle & \dots & \langle y_2 x_N \rangle & C_{2N}^{YY} & \langle y_2 z_N \rangle \\ \dots & \dots & \dots & \dots & \dots & \dots & \langle x_N^2 \rangle & \langle x_N y_N \rangle & \langle x_N z_N \rangle \\ \dots & \dots & \dots & \dots & \dots & \dots & \langle y_N^2 \rangle & \langle y_N y_N \rangle & \langle y_N z_N \rangle \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \langle x_N^2 \rangle - \langle y_N^2 \rangle & 1 - \langle x_N^2 \rangle - \langle y_N^2 \rangle \end{pmatrix} \ge 0 \quad (7) $$

Notice that we mark in blue the entries replaced with the available data. All the other entries (black terms $\langle \cdots \rangle$) are unknowns, which represent unobserved correlations over the $p$ distribution. If other correlators were known (e.g., $C_{12}^{XY}$), they would simply replace the corresponding free variables in Eq. (7) (namely $\langle x_1 y_2 \rangle$): reducing the number of free variables makes it harder to complete the matrix $\Gamma \ge 0$, and therefore makes it easier to detect entanglement. Notice that correlations such as $\langle x_i y_i \rangle$ have no experimental meaning in quantum physics, as they involve the simultaneous measurement of two incompatible observables, namely $\hat{X}_i$ and $\hat{Y}_i$ on the same qubit. However, they represent perfectly well-defined 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 $\Gamma$ matrix with such unobserved correlations, such that $\Gamma \ge 0$. Crucially, solving this problem is a so-called *semidefinite program* [25], for which efficient convex-optimization algorithms are available. As we illustrate in Sec. 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 $\Gamma \ge 0$ cannot be completed, the theory of semidefinite programming allows one to derive an entanglement witness of the form:

$$ \sum_{r=1}^{R} w_r C_r \stackrel{\hat{\rho}_{\text{sep}}}{\le} 1 - \lambda, \quad (8) $$

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

$$ \sum_{r=1}^{R} w_r C_r = 1, \quad (9) $$

ultimately proving that the quantum state generating these data is entangled. The parameter $\lambda > 0$ in Eq. (8) can be interpreted as the *noise robustness* of the data. Indeed, if the quantum state $\hat{\rho}$ generating the data is mixed with white noise: $\hat{\rho} \rightarrow (1-\lambda)\hat{\rho} + \lambda \mathbb{1}/D$ with $D = 2^N$ the dimension of the Hilbert space, then using the fact that Pauli observables are traceless, we have $\{C_i^a, C_{ij}^{ab}\} \rightarrow \{(1-\lambda)C_i^a, (1-\lambda)C_{ij}^{ab}\}$. The parameter $\lambda$ thus exactly quantifies the maximal amount of white noise, which can be tolerated before entanglement detection becomes impossible with Eq. (7). Hence, by using a SDP to minimize the noise strength $\lambda$ for which the noisy data becomes compatible with a $\Gamma \ge 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 variables $\mathbf{m}' = \{1\} \cup \{a_i\} \cup \{a_i b_j\} \cup \{a_i b_j c_k\} \cup \dots$, where $a_i, b_j, c_k$ are any components of the local classical variables $\{\mathbf{n}_i\}$. One then constructs the (PSD) correlation matrix $\Gamma'_{\alpha,\beta} = \langle m'_{\alpha} m'_{\beta} \rangle$ over the $p$ distribution. Verifying the compatibility of the data with $\Gamma' \ge 0$ is again a semidefinite program, which can be solved at a computational (memory) cost scaling at most as $O[\text{length}(\mathbf{m}')^2]$. The matrix $\Gamma$ [Eq. (7)] is obtained as a submatrix of $\Gamma'$, and therefore the condition $\Gamma' \ge 0$ is stronger than $\Gamma \ge 0$. Including in $\mathbf{m}'$ all monomials up to degree $l = 1, 2, 3, \dots$, one defines a systematic hierarchy of positive-semidefinite conditions, which are necessarily obeyed if the underlying state $\hat{\rho}$ 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 $\Gamma' \ge 0$—this property is a consequence of the variables $\{\mathbf{n}_i\}$ being compact: $|\mathbf{n}_i| = 1$ for all $i$. Therefore, the approach presented here defines, in the limit $l \to \infty$, 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)^2]$ 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[\{\mathbf{n}_i\}]$ [25,38]. Notice that in practice, the computational cost of higher-level tests ($l \ge 2$) increases rapidly, especially for hundreds of qubits. However, we provide in Sec. 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 PT criteria [33,34,39,40]: a state $\hat{\rho}$ is compatible with Eq. (7) if and only if the state $\hat{\rho}^{\text{PT}}$ is compatible with Eq. (7), where $\hat{\rho}^{\text{PT}}$ is obtained by applying partial transposition on any subset of qubits. Indeed, PT leaves invariant the Pauli matrices $\hat{X}_i$ and $\hat{Z}_i$, while $\hat{Y}_i$ is changed into $-\hat{Y}_i$. Therefore, all correlations involving $\hat{Y}_i$ are changed into their opposite. The corresponding matrix $\Gamma^{\text{PT}}$ is then obtained from $\Gamma$ [Eq. (7)] by a simple change of basis, in which $(x_i, y_i, z_i)$ is changed into $(x_i, -y_i, z_i)$ for the qubits where PT is applied. Therefore, $\Gamma^{\text{PT}}$ can be completed as a PSD matrix if and only if $\Gamma$ can be completed as a PSD matrix. As discussed in Sec. 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 [41], namely a statistical mixture of white noise with a spin singlet:

$$ \hat{\rho}_{\lambda} = (1 - \lambda)|\Psi\rangle\langle\Psi| + \frac{\lambda}{4}1, \quad (10a) $$

$$ |\Psi\rangle = \frac{1}{\sqrt{2}}(|\uparrow\downarrow\rangle - |\downarrow\uparrow\rangle), \quad (10b) $$

where $0 \le \lambda \le 1$. The state $\hat{\rho}_{\lambda}$ is separable if and only if $\lambda \ge 2/3$. Let us show that Eq. (7) is tight for the Werner state, namely that it detects entanglement whenever $\lambda < 2/3$. The Werner state is SU(2) invariant, and one finds $C_1^a = C_2^a = 0$ (for $a \in \{X, Y, Z\}$) and $C_{12}^{ab} = -\delta_{ab}(1 - \lambda)$. However, this detailed property, impossible 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_{12}^{XX} + C_{12}^{YY} + C_{12}^{ZZ}$ 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

$$ \Gamma = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ . & 1/3 & 0 & 0 & \mathbf{c}/3 & 0 & 0 \\ . & . & 1/3 & 0 & 0 & \mathbf{c}/3 & 0 \\ . & . & . & 1/3 & 0 & 0 & \mathbf{c}/3 \\ . & . & . & . & 1/3 & 0 & 0 \\ . & . & . & . & . & 1/3 & 0 \\ . & . & . & . & . & . & 1/3 \end{pmatrix} \ge 0. \quad (11) $$

Reorganizing the lines and columns, the matrix $\Gamma$ is block diagonal, and is PSD if and only if all blocks are PSD, that is, if and only if $(\frac{1}{c} \begin{smallmatrix} c \\ 1 \end{smallmatrix}) \ge 0$, if and only if $|c| \le 1$. This establishes that entanglement is detected whenever $|C_{12}^{XX} + C_{12}^{YY} + C_{12}^{ZZ}| > 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 $|\uparrow\uparrow\rangle$ is such that $C_{12}^{XX} + C_{12}^{YY} + C_{12}^{ZZ} = C_{12}^{ZZ} = 1$. In the case of the Werner state, we have $c = -3(1 - \lambda)$, from which we recover the known result that the state, Eq. (10), is entangled for $\lambda < 2/3$. Notice that Eq. (7) is not tight for all two-qubit states: by producing random two-qubit states, we could find entangled states (as detected by the concurrence criterion [42]), which are nevertheless compatible with a PSD correlation matrix as in Eq. (7).

### IV. ROBUST DETECTION OF ENTANGLEMENT IN MANY-BODY SYSTEMS

We then choose to benchmark our entanglement detection method on paradigmatic lattice quantum spin models. Motivated by an ultracold atoms' experiment [43] realized a few years ago [44], we first focus on the entanglement generated by a single impurity propagating along a one-dimensional XX chain (Sec. IVA). As a main result, we find 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. [44]. We then consider a two-dimensional system, where entanglement is generated by a XX Hamiltonian with $1/r^3$ interactions, from an initial state with all spin polarized along X (Sec. IV B). This example of a many-body quench with power-law interactions is especially motivated by Rydberg arrays [18,20,21,45–49], ultracold magnetic atoms [50–52], nitrogen-vacancy centers in diamond [53,54], and trapped-ion systems [11,13,55–60], where related spin Hamiltonians have been implemented. In this case, our algorithm leads 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 explore 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 (Sec. 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 $\rho$ with white noise, modeled as a completely mixed state: $\hat{\rho} \rightarrow (1 - \lambda)\hat{\rho} + \lambda \mathbb{1}/2^N$. The noise robustness $\lambda^*$ is then defined as the value of $\lambda$ 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 $\lambda^*$ (cf. Sec. II and Appendix A2 for details). To generate the numerical SDP problems, we use the software Ncpol2sdpa [61], and we solve the SDPs with Mosek [62]. We release an open source code [63], which allows one to recover the results of Sec. IVA, 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-neighbor interactions:

$$ \hat{\mathcal{H}}_{XX} = -J \sum_{i=0}^{N-1} [\hat{X}_i \hat{X}_{i+1} + \hat{Y}_i \hat{Y}_{i+1}], \quad (12) $$

with $J = 1$ a global energy scale, and with periodic boundary conditions. As initial state, we consider the ferromagnetic state $|\Psi_0\rangle = \otimes_i |\uparrow\rangle$. We assume that at time $t = 0$, the spin at $i = 0$ is flipped into $|\downarrow\rangle$ [44,64]. This central excitation then propagates along the chain under the XX Hamiltonian. As $\hat{\mathcal{H}}_{XX}$ conserves the total magnetization along Z, the dynamics occurs in the N-dimensional manifold of states generated by $\{\hat{\sigma}_i^- | \Psi_0\rangle\}_{i=0}^{N-1}$ [with $\hat{\sigma}_i^- = (\hat{X}_i - i\hat{Y}_i)/2$ 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. [44], the propagation of entanglement was observed through a lower bound to the pairwise concurrence [64], which measures the entanglement of the two-body reduced state $\hat{\rho}_{ij}$ [42]. Here, our main result is that using the same information as in the experiment of Ref. [44] [namely, the transverse correlations $C_{ij}^{\perp} := (C_{ij}^{XX} + C_{ij}^{YY})/2$, the magnetization $C_{ij}^{Z}$ and the longitudinal correlations $C_{ij}^{ZZ}$], 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

$$ \varphi_r(t) = N^{-1} \sum_{k=0}^{N-1} \exp \left[ \frac{2i\pi kr}{N} + it \cos \left( \frac{2\pi k}{N} \right) \right], \quad (13a) $$

$$ C_i^Z = 1 - 2|\varphi_i|^2, \quad (13b) $$

$$ C_{ij}^{\perp} = 2\mathrm{Re}(\varphi_i^* \varphi_j), \quad (13c) $$

$$ C_{ij}^{ZZ} = 1 - 2(|\varphi_i|^2 + |\varphi_j|^2). \quad (13d) $$

As discussed in Appendix A1, in order to implement our algorithm, we may use the symmetries of the problem to drastically reduce the number of nonzero variables in Eq. (7), greatly improving the scalability. The resulting witness at each time, reconstructed via the algorithm described in Appendix A2, 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 $\hat{W} = \sum_i w_i^Z \hat{Z}_i + \sum_{i \neq j} [w_{ij}^{ZZ} \hat{Z}_i \hat{Z}_j + w_{ij}^{\perp} (\hat{X}_i \hat{X}_j + \hat{Y}_i \hat{Y}_j)/2]$. In 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 choose 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 [64], as measured in the experiment of Ref. [44]. The SDP witness is about one order of magnitude more robust against 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}^{\perp}$ [Fig. 2(e)], the qubits whose contribution to the witness is the largest are located close to $i = -j = \pm 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
<figure><img src="image_4.png" /><figcaption>FIG. 2. Single spin flip in a XX chain with $N = 64$ spins (cf. experiment of Ref. [44]). (a)–(c) One- and two-body correlations at time $tJ = 10$, used as input to our SDP algorithm [for the sake of readability, on (c) we plot the connected correlations $C_{ij}^{ZZ} - C_i^Z C_j^Z$]. (d)–(f) Corresponding coefficients in the optimal entanglement witness. (g) On a semilogarithmic scale, noise robustness of the violation of the witness as a function of time (blue dots); and theoretical prediction for the noise robustness of the concurrence lower-bound measured in Ref. [44] (solid orange line).</figcaption></figure>

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:

$$ \hat{\mathcal{H}}_{XX} = J \sum_{1 \le i < j \le N} \frac{\hat{X}_i \hat{X}_j + \hat{Y}_i \hat{Y}_j}{r_{ij}^3} + h \sum_{i=1}^{N} \hat{X}_i, \quad (14) $$

where $r_{ij}$ denotes the distance between spins $i$ and $j$, arranged over a $N = L \times 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 [49], and to trapped ions [60]. As initial state, we consider a ferromagnetic state along $X$: $|\Psi_0\rangle = \otimes_i |+\rangle_i$ with $|+\rangle = (|\uparrow\rangle + |\downarrow\rangle)/\sqrt{2}$. For this particular initial state, the dynamics is invariant under the change $\hat{\mathcal{H}}_{XX} \to -\hat{\mathcal{H}}_{XX}$; and

$|\Psi_0\rangle$ represents the mean-field ground state of $-\hat{\mathcal{H}}_{XX}$ [65]. The dynamics is then well approximated by a semiclassical spin-wave approach [66], involving bosonic Gaussian states, whose stability is further enhanced by introducing the symmetry-breaking term $h \sum_i \hat{X}_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 select 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 use the one-body terms $C_i^X$ (for all qubits $i$) and the two-body terms $C_{ij}^{aa}$ (for $a = X, Y, Z$ and all pairs $i < j$). Once again, we use symmetries to reduce the number of free variables in the SDP algorithm (see Appendices A1, and A2 for details on the algorithm used to reconstruct an entanglement witness from the data). We first consider systems with periodic boundary conditions, such that the data are translationally invariant (TI). In this case, we can analyze analytically the witnesses found by our algorithm, and
generalize them to a whole family of entanglement witnesses in its own right. We then consider 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 find the following family of witnesses:

$$ \sum_{a \in \{X,Y,Z\}} \sum_{j \neq j'} e^{i[\phi_a(j') - \phi_a(j)]} C_{jj'}^{aa} \geq -N, \quad (15) $$

where $\phi_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 (Sec. II). Using that $N = \sum_i (x_i^2 + y_i^2 + z_i^2)$, we then have $\sum_{a \in \{X,Y,Z\}} \sum_{j \neq j'} e^{i[\phi_a(j') - \phi_a(j)]} C_{jj'}^{aa} + N^{\text{sep}} = \sum_{a \in \{X,Y,Z\}} \langle |\sum_j e^{-i\phi_a(j)} a_j|^2 \rangle \geq 0$.

The witness of Eq. (15) turns out to be very similar to existing results in the literature [67–69]. There are however two important differences: on the technical side, the witness of Eq. (15) involves local phases $\phi_a(j)$, which may depend on the spin direction $a$ (this possibility was not pointed out in the mentioned Refs. [67–69]); 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 $\phi_a(j) = k_a \cdot r_j$ with $r_j$ the position of the $j$th subsystem. This leads to the structure factor $S_k^a = N^{-1} \sum_{j,j'} e^{i(k \cdot (r_j - r_{j'})} C_{jj'}^{aa}$ (notice that we include the $j = j'$ term in the summation, corresponding to a term $C_{jj'}^{aa} = 1$). Although we discover 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

$$ S_{k_X}^X + S_{k_Y}^Y + S_{k_Z}^Z \geq 2. \quad (16) $$

Notice that we define 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 wave vector $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

$$ W_{\text{opt}} = \min_{\mathbf{k}} S_{\mathbf{k}}^{X} + \min_{\mathbf{k}} S_{\mathbf{k}}^{Y} + \min_{\mathbf{k}} S_{\mathbf{k}}^{Z} \geq 2. \quad (17) $$

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 then generalize this result to Eq. (15). The witness Eq. (17) is then violated if the fluctuations at these optimal wave vectors 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 wave vectors may vary over time in different ways for different spin components; and therefore the witnesses detect entanglement, independently of the specific SDP algorithm we use 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_k^a = N^{-1} \sum_{j,j'} e^{i(k \cdot (r_j - r_{j'})} \text{Tr}[\hat{S}_j^a \hat{S}_{j'}^a]$ with $\hat{S}_j^a$ the spin-s observable in direction $a$ for subsystem $j$. Using that $Ns^2 \geq \sum_{i=1}^N \sum_{a \in \{X,Y,Z\}} |\text{Tr}[\hat{S}_j^a]|^2$, one easily shows that $S_{k_X}^X + S_{k_Y}^Y + S_{k_Z}^Z \geq Ns$ for all separable states. This generalizes a result of Ref. [70] to arbitrary components of the structure factors.

**Numerical results.** In Fig. 3(a), we plot the structure factor at time $tJ = 1$, for wave vectors $k = (k, k)$, in a $20 \times 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. As we consider open boundary conditions, the correlations have no translation invariance, and as a consequence the resulting witness loses this symmetry too. For comparison, we also plot the noise robustness of the TI witness of Eq. (17), 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 cannot find an analytical expression for these non-TI data-driven witnesses. By a variational search over separable states, we can, however, verify that the separable bound obtained by our SDP algorithm is 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_i^a$ and only inter-AB correlations $C_{ij}^{ab}$ where $i \in A$ and $j \in B$.
<figure><img src="image_5.png" /><figcaption>FIG. 3. Quench in the two-dimensional XX model on a 20 × 20 square lattice with open boundary conditions [Eq.(14)]. (a) Structure factor at time t = 1/J. The circles mark the wave vectors at which the structure factor for the X, Y, and Z components is minimal. This leads to the TI witness of Eq. (17). (b) Blue dots: noise robustness of the optimal witness found by our SDP algorithm using as input the magnetization C_i^X and two-body correlations C_{ij}^{XX}, C_{ij}^{YY}, and C_{ij}^{ZZ}, with an initial state polarized along x. Orange line: noise robustness of the best TI witness of Eq. (17), whose reconstruction is illustrated on (a) for time t = 1/J.</figcaption></figure>

It is straightforward to notice that if non-full-separability can be proved from this knowledge, then the state must be bipartite entangled. Indeed, if the AB state is bipartite separable: $\hat{\rho}_{AB}^{\text{bisep}} = \sum_k p_k \hat{\rho}_A^{(k)} \otimes \hat{\rho}_B^{(k)}$, then we may define $\hat{\rho}_{AB}^{\text{fullsep}} = \sum_k p_k \otimes_{i=1}^N \hat{\rho}_i^{(k)}$ with $\hat{\rho}_i^{(k)} = \text{Tr}_{j \neq i}[\hat{\rho}_A^{(k)} \otimes \hat{\rho}_B^{(k)}]$. One can verify that if $\hat{\rho}_{AB}^{\text{bisep}}$ reproduces the data, so does $\hat{\rho}_{AB}^{\text{fullsep}}$; 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 Heisenberg chain is described by the Hamiltonian:

$$ \hat{\mathcal{H}}_{\text{Heis}} = (J/4) \sum_{i=1}^{N} [\hat{X}_i \hat{X}_{i+1} + \hat{Y}_i \hat{Y}_{i+1} + \hat{Z}_i \hat{Z}_{i+1}], \quad (18) $$

and transverse-field Ising chain by

$$ \hat{\mathcal{H}}_{\text{Ising}} = -(J/4) \sum_{i=1}^{N} [\hat{Z}_i \hat{Z}_{i+1} + g \hat{X}_i]. \quad (19) $$

In Eqs. (18) and (19), we assume periodic boundary conditions. The parameter J is an overall energy scale (set to J = 1 in our computations), and g in Eq. (19) is the transverse-field amplitude. We consider thermal states $\hat{\rho} = Z^{-1} \exp[-\hat{\mathcal{H}}/T]$ with T the temperature and $Z = \text{Tr}(\exp[-\hat{\mathcal{H}}/T])$ the partition function. As input to

our algorithm, we use all one- and two-body correlations, which are invariant under translations. These data were computed with quantum Monte Carlo. The Heisenberg model being SU(2) invariant, one-body terms vanish, and two body-terms are of the form $C_{ij}^{ab} = C_{|i-j|} \delta_{ab}$ (with a, b ∈ {X, Y, Z}). For the Ising model the symmetries imply $C_i^Y = C_i^Z = 0$, and $C_{ij}^{ab} = C_{|i-j|}^a \delta_{ab}$ (namely, off-diagonal correlations $a \neq b$ vanish).

**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 notice that entanglement is detected if and only if the nearest-neighbor two-body reduced density matrix $\rho_{N/2-1,N/2}$ is itself entangled (as detected by the concurrence [42]). While illustrating the relatively short-range nature of entanglement in these thermal states, we cannot go beyond the mere witnessing of entanglement among nearest neighbors. We therefore consider a partitioning maximizing the AB interface, that is A|B|A|B|A|B|... (namely, A = {0, 2, 4, ..., N − 2} and B = {1, 3, 5, ..., N − 1}). This leads us to discover new bipartite entanglement witnesses. Similarly to the case of TI multipartite entanglement witnesses [see Eq. (15)], we can analytically characterize them, and extend them to a full family of witnesses. We define

$$ W_a = \sum_{i \in A} \sum_{j \in B} K_{j-i} C_{ij}^{aa} \cos[\phi_a(i) - \phi_a(j)], \quad (20) $$

where $\phi_a(j)$ are arbitrary local phases, and with coefficients given by

$$ K_r = \frac{2}{N} \sum_{k=-\frac{N}{4}+1}^{\frac{N}{4}-1} \exp\left(\frac{2i\pi}{N}kr\right) \quad (21) $$

$$ = \frac{2}{N} \left[ \frac{\sin(\pi r/2)}{\tan(\pi r/N)} - \cos\left(\frac{\pi r}{2}\right) \right] \quad (22) $$

$$ \stackrel{N \gg r}{\sim} \frac{2}{\pi r} \sin\left(\frac{\pi r}{2}\right). \quad (23) $$

Notice that for $i \in A$ and $j \in 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(\pi r/N)$. The violated witnesses are then of the form

$$ W = W_X + W_Y + W_Z \geq -\frac{N}{2}. \quad (24) $$

The proof of the witness inequality (24) is given in Appendix B. For the Heisenberg chain, the optimal choice of the phases is $\phi_a(i) = 0$. For the Ising chain, it is $\phi_a(i) = 0$ for $i \in A$, and $\phi_X(j) = \phi_Z(j) = \pi$ for $j \in B$, and $\phi_Y(j) = 0$.

As illustrated in Fig. 4, the bipartite entanglement witnesses of Eq. (24) allow one to detect entanglement in
<figure><img src="image_6.png" /><figcaption>FIG. 4. Entanglement in (a) the Heisenberg chain [Eq. (18)] and (b) the quantum Ising chain at $g = 1$ [Eq. (19)], for $N = 64$ spins. The bipartition is of the form $A|B|A|B|A|B...$. Solid orange line: (minus) the concurrence [42] between nearest neighbors. Blue dots: bipartite witness $1 + (2W/N)$ with $W$ given in Eq. (24). Dashed-dotted green line: multipartite witness $S_{k_X}^X + S_{k_Y}^Y + S_{k_Z}^Z - 2$ at the optimal wave vectors $k_a$ [see text and Eq. (17)] (for all witnesses, entanglement is detected for negative values).</figcaption></figure>

regimes where all two-body reduced density matrices are separable (as measured by the concurrence [42]).

**Comparison with the structure factor witnesses.** For the sake of completeness, we also evaluate the multipartite entanglement witness based on the structure factor [Eq. (17)] discussed in Sec. 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 wave vectors $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 = \pi$ in the structure factor. Concomitantly, fluctuations of the uniform magnetization (at $k = 0$) are suppressed, reaching a spin-singlet state in the ground state ($S_0^a = 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 $\sum_a S_0^a \ge 2$. As previously reported in Ref. [36], it is violated up to a temperature $T/J \approx 1.4$ [Fig. 4(a)]. In contrast, at low temperature and around the quantum critical point $g = 1$ [71], the Ising model develops ferromagnetic correlations for the *Z* component of the spin, and the structure factor is minimal at $k_Y = \pi$. The uniform magnetization along *Y* is slightly squeezed below the standard quantum limit [72], and the structure factor is minimal at $k_Y = 0$. Finally, correlations in the *X* direction (the direction of the transverse field) are strongly ferromagnetic, and are suppressed at $k_X = \pi$. We find 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_\pi^X + S_0^Y + S_\pi^Z \ge 2$

is violated for temperatures $T/J \lesssim 0.38$. In contrast a criterion based on the quantum Fisher information (based on the dynamical structure factor for ZZ correlations [73], and which is considerably more challenging to estimate, both in theory and in experiments [74,75]) is violated for $T/J \lesssim 0.11$; and a witness based on the same data as used in the present work [36], and optimized for $T/J = 0.28$, is violated for $T/J \lesssim 0.31$. The entanglement 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., Eqs. (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 [27] 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 [76].

Several research directions are now open to future works. First, 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 [77], where entanglement could be revealed by string-order-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 [15–21], and where entanglement is commonly assumed to be an essential resource [2].

## ACKNOWLEDGMENTS

We thank Tommaso Roscilde for providing us the quantum Monte Carlo data used in Section IV.C. We thank Felix Huber for drawing our attention to Ref. [24] 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 No. 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).

## APPENDIX A: TECHNICAL CONSIDERATIONS

### 1. Symmetrizing the unknown $p$ distribution

In this section, we explain how the number of free variables in the SDP algorithm solving Eq. (7) can be drastically reduced by symmetrizing the underlying $p(\{x_i, y_i, z_i\})$ distribution. This symmetrization relies only on the nature of the data themselves, and does not assume that the quantum state has specific symmetries.

(1) First, let us assume that the available data consist of $\{C_i^Z, C_{ij}^{XX}, C_{ij}^{YY}, C_{ij}^{ZZ}\}$, namely the magnetization along the direction Z for all qubits, and pair correlations along the

same direction. The data we used in analyzing the quench in the two-dimensional XX model with $1/r^3$ interactions, and thermal states of the quantum Ising model, have this structure (exchanging the role of the Z and X direction). Let us assume that a distribution $p_0$ exists, which reproduces the data. One may then consider the distributions $p_1 = p_0(\{-x_i, y_i, z_i\})$, $p_2 = p_0(\{x_i, -y_i, z_i\})$, and $p_3 = p_0(\{-x_i, -y_i, z_i\})$. Clearly, one has that $\langle z_i \rangle_{p_0} = \langle z_i \rangle_{p_k}$ for $k \in \{1, 2, 3\}$, and similarly for two-body terms: $\langle a_i a_j \rangle_{p_0} = \langle a_i a_j \rangle_{p_k}$ for $a \in \{x, y, z\}$ and $k \in \{1, 2, 3\}$. In other words, $p_0$ reproduces the data if and only if $p_k$ reproduces the data for $k \in \{1, 2, 3\}$. One may then consider the distribution $p_4 = (p_0 + p_1 + p_2 + p_3)/4$, which also reproduces the data. Importantly, $p_4$ is such that $\langle x_i \rangle_{p_4} = 0$ (since under $p_4$, $x_i$ has the same probability as $-x_i$). Similarly, $\langle y_i \rangle_{p_4} = 0$, and two-body terms are such that $\langle x_i y_j \rangle_{p_4} = \langle x_i y_j \rangle_{p_4} = \langle y_i z_j \rangle_{p_4} = 0$. Therefore, without loss of generality, one may impose in the SDP algorithm that $\langle x_i \rangle = \langle y_i \rangle = \langle x_i y_j \rangle = \langle x_i z_j \rangle = \langle y_i z_j \rangle = 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 $\langle x_i^2 \rangle$ and $\langle y_i^2 \rangle$ for $i \in [N]$).

(2) Let us then consider a situation where the data consist of $\{C_i^Z, C_{ij}^{XX} + C_{ij}^{YY}, C_{ij}^{ZZ}\}$, 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 $\langle x_i x_j + y_i y_j \rangle_p = \langle x_i x_j + y_i y_j \rangle_{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 $\langle x_i x_j \rangle_q = \langle y_i y_j \rangle_q$ for all pairs $(i, j)$, we may assume $C_{ij}^{XX} = C_{ij}^{YY}$ 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 $\langle x_i x_j \rangle_p = \langle x_i x_j \rangle_{p''}$. Considering $(p + p'')/2$, we may impose $\langle x_i \rangle = 0$, and also $\langle x_i y_j \rangle = \langle x_i z_j \rangle = 0$ for all pairs $(i, j)$. The same reasoning apply interchanging the role of x and y, leading to $\langle y_i \rangle = \langle y_i z_j \rangle = 0$. In conclusion, we may solve the SDP in the form of Eq. (7) imposing that all terms are zero, except the data $C_i^Z, C_{ij}^{XX} = C_{ij}^{YY}$ [obtained as $(C_{ij}^{XX} + C_{ij}^{YY})/2$], and $C_{ij}^{ZZ}$; and except the diagonal terms with $\langle x_i^2 \rangle = \langle y_i^2 \rangle$. This reduces the number of free variables from order $O(N^2)$ to $N$ (for instance, all $\langle z_i^2 \rangle$ for $i \in [N]$, greatly improving the scalability).

(3) Finally, if one only knows $c_{ij} := C_{ij}^{XX} + C_{ij}^{YY} + C_{ij}^{ZZ}$, by a similar reasoning one may consider a distribution $p(\{x_i, y_i, z_i\})$, which is invariant under all rotations, satisfying $\langle x_i \rangle = \langle y_i \rangle = \langle z_i \rangle = 0$, $\langle x_i y_j \rangle = \langle x_i z_j \rangle = \langle y_i z_j \rangle = 0$, and $\langle x_i x_j \rangle = \langle y_i y_j \rangle = \langle z_i z_j \rangle$. In this case, there is no free variable at all in Eq. (7) (all two-body diagonal terms are $1/3$, and the only nonzero off-diagonal entries are $C_{ij}^{XX} = C_{ij}^{YY} = C_{ij}^{ZZ}$, 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 \neq 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.

## 2. Robust entanglement witness from one- and two-body correlations on N qubits

Here we provide details on the way to extract an entanglement witness in the form of Eq. (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_i^\alpha$ and two-body terms $C_{ij}^{ab}$ are measured. We denote generically these data as $\{C_\alpha\}_{\alpha=1}^R := \{C_i^\alpha, C_{ij}^{ab}\}$, where $\alpha$ 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 $\lambda$ such that $\{(1 - \lambda)C_\alpha\}$ is compatible with a PSD correlation matrix, as explained in the main text [Sec. II and, in particular, Eq. (7)]. That is, we aim at solving the problem:

$$ \min_{\lambda \ge 0} \lambda \quad \text{such that,} \\ \text{(PSD)} \quad \Gamma \ge 0, \qquad \text{(A1a)} $$

$$ \text{(data)} \quad \Gamma_{(n,m)(\alpha)} = (1 - \lambda)C_{\alpha} \quad (\alpha \in [R]), \qquad \text{(A1b)} $$

$$ \text{(Pauli)} \quad \sum_{a=1}^{3} \Gamma_{3(i-1)+a,3(i-1)+a} = 1 \quad (i \in [N]). \qquad \text{(A1c)} $$

The (symmetric) matrix $\Gamma = (\Gamma_{n,m})_{0 \le n,m \le 3N}$ is the correlation matrix for the variables $\{1\} \cup \{x_i, y_i, z_i\}_{i=1}^N$ parametrizing separable states of $N$ qubits (Sec. II), and is therefore PSD [Eq. (A1a)]. In Eq. (A1b), we constrain the relevant entries of the $\Gamma$ matrix to reproduce the data (with a $1 - \lambda$ noise prefactor); we introduce $(n, m)(\alpha)$ to denote the pair of indices $(n, m)$ containing the data $C_\alpha$. Specifically [see Eq. (7)], $C_i^\alpha$ is contained in $\Gamma_{0,3(i-1)+a}$, and $C_{ij}^{ab}$ in $\Gamma_{3(i-1)+a,3(j-1)+b}$. Finally, Eq. (A1c) enforces $x_i^2 + y_i^2 + z_i^2 = 1$ in the $\Gamma$ matrix; as this condition descends from properties of the Pauli matrices (Sec. II), we call it the Pauli constraint.

**Standard primal form.** We now rewrite in its so-called standard primal form [25] the SDP defined in Eq. (A1):

$$ \begin{align} & \min_{X \ge 0} \langle M, X \rangle \quad \text{such that,} \tag{data} \\ & \text{(data)} \quad \langle A_{\alpha}^{\text{data}}, X \rangle = C_{\alpha} \quad (\alpha \in [R]), \nonumber \\ & \text{(Pauli)} \quad \langle A_{i}^{\text{Pauli}}, X \rangle = 1 \quad (i \in [N]). \tag{A2} \end{align} $$

We introduce the matrix scalar product $\langle X, Y \rangle = \text{Tr}(X^T Y) = \sum_{ij} X_{ij} Y_{ij}$. The involved matrices have the following block-diagonal form:

$$ M = \left( \begin{array}{c|c} 1 & 0 \\ \hline 0 & 0 \end{array} \right), X = \left( \begin{array}{c|c} \lambda & 0 \\ \hline 0 & \Gamma \end{array} \right), \qquad \text{(A3a)} $$

$$ A_{\alpha}^{\text{data}} = \left( \begin{array}{c|c} C_{\alpha} & 0 \\ \hline 0 & \Sigma_{\alpha}^{\text{data}} \end{array} \right), A_{i}^{\text{Pauli}} = \left( \begin{array}{c|c} 0 & 0 \\ \hline 0 & \Sigma_{i}^{\text{Pauli}} \end{array} \right), \qquad \text{(A3b)} $$

where the $\Sigma$ matrices are given by

$$ (\Sigma_{\alpha}^{\text{data}})_{n,m} = \begin{cases} 1 & \text{if } (n,m) = (n,m)(\alpha) \\ 0 & \text{otherwise} \end{cases}, \qquad \text{(A4a)} $$

$$ (\Sigma_i^{\text{Pauli}})_{n,m} = \sum_{a=1}^{3} \delta_{n,m} \delta_{m,3(i-1)+a}. \qquad \text{(A4b)} $$

The SDP, Eq. (A2), is equivalent to the problem, Eq. (A1), of finding the minimal noise $\lambda$ for which the $\Gamma$ matrix in Eq. (7) becomes positive semidefinite, for a set of noisy data $\{(1 - \lambda)C_i^\alpha, (1 - \lambda)C_{ij}^{ab}\}$. Therefore, an optimal solution $\lambda^* > 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 [25] corresponding to Eq. (A2):

$$ \begin{gathered} \max_{\vec{w}, \vec{w}^{\text{Pauli}}} \sum_{\alpha=1}^{R} w_{\alpha} C_{\alpha} + \sum_{i=1}^{N} w_{i}^{\text{Pauli}} \quad \text{such that,} \\ \sum_{\alpha=1}^{R} w_{\alpha} A_{\alpha}^{\text{data}} + \sum_{i=1}^{N} w_{i}^{\text{Pauli}} A_{i}^{\text{Pauli}} \le M. \end{gathered} \qquad \text{(A5)} $$

Using the expressions of the matrices $M$ and $A$'s [Eqs. (A3) and (A4)], we rewrite Eq. (A5) as

$$ \begin{align} & \max_{\vec{w}, \vec{w}^{\text{Pauli}}} \vec{w} \cdot \vec{C} + \sum_{i=1}^{N} w_i^{\text{Pauli}} \quad \text{such that,} \nonumber \\ & \vec{w} \cdot \vec{C} \le 1 \nonumber \\ & \sum_{\alpha=1}^{R} w_{\alpha} \Sigma_{\alpha}^{\text{data}} + \sum_{i=1}^{N} w_i^{\text{Pauli}} \Sigma_i^{\text{Pauli}} \le 0. \tag{A6} \end{align} $$

We denote $\{\vec{w}^*, (\vec{w}^{\text{Pauli}})^*\}$ the optimal solution to this dual problem.

**Strong duality.** As for all semidefinite programs [25], any primal feasible $X$ [that is, a PSD matrix $X \ge 0$ satisfying the constraints in Eq. (A2)] yields an upper bound
to any dual feasible [that is, any {$w_\alpha, w_i^{\text{Pauli}}$} satisfying the constraints in Eq. (A5)]. Indeed, for any such feasible $X$ and $\vec{w}$, we have $\langle M, X \rangle - \sum_\alpha w_\alpha C_\alpha - \sum_i w_i^{\text{Pauli}} = \langle (M - \sum_\alpha w_\alpha A_\alpha^{\text{data}} - \sum_i w_i^{\text{Pauli}} A_i^{\text{Pauli}}), X \rangle$, where we used the constraints in Eq. (A2). Using then the constraint in Eq. (A5), we notice that $(M - \sum_\alpha w_\alpha A_\alpha^{\text{data}} - \sum_i w_i^{\text{Pauli}} A_i^{\text{Pauli}}) \ge 0$. As $X \ge 0$, and using the fact that $\langle Y, X \rangle \ge 0$ for any two PSD matrices $Y$ and $X$, we conclude that $\langle M, X \rangle - \sum_\alpha w_\alpha C_\alpha - \sum_i w_i^{\text{Pauli}} \ge 0$. In particular, the primal optimum upper bounds the dual optimum: $\lambda^* \ge \vec{w}^* \cdot \vec{C} + \sum_i (w_i^{\text{Pauli}})^*$, a property known as weak duality [25]. 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 [25] is that both primal and dual problems are strictly feasible. The primal problem is strictly 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 $\lambda = 1$ and $\Gamma = 1/3$ in Eq. (A3) [or, equivalently, in Eq. (A1)]. The dual problem is strictly feasible if one can exhibit some {$w_\alpha, w_i^{\text{Pauli}}$} satisfying the constraints in Eq. (A6) as strict inequalities: such strictly feasible $w$'s are readily obtained as $w_\alpha = 0$ and $w_i^{\text{Pauli}} = -1$ (such that $\sum_i w_i^{\text{Pauli}} \Sigma_i^{\text{Pauli}} = -1$). Two important properties follow from strong duality: (1) $\lambda^* = \vec{w}^* \cdot \vec{C} + \sum_i (w_i^{\text{Pauli}})^*$; (2) as a consequence, we also have $\langle [M - \sum_\alpha w_\alpha^* A_\alpha^{\text{data}} - \sum_i (w_i^{\text{Pauli}})^* A_i^{\text{Pauli}}], X^* \rangle = 0$. Using the expression of $M$ and the $A$'s matrices [Eq. (A3)], this implies, in particular, $(1 - \vec{w}^* \cdot \vec{C})\lambda^* = 0$. Whenever entanglement is detected ($\lambda^* > 0$), we have therefore $\vec{w}^* \cdot \vec{C} = 1$, so that $-\sum_i (w_i^{\text{Pauli}})^* = 1 - \lambda^*$.

**Entanglement witness.** The coefficients $w_\alpha^*$ define an entanglement witness whose separable bound is given by $1 - \lambda^*$:

$$ \sum_{\alpha} w_{\alpha}^{*} C_{\alpha}^{\text{sep}} \leq 1 - \lambda^{*}, \qquad (A7) $$

which holds for all separable data {$C_\alpha^{\text{sep}}$} (namely, data compatible with a separable state), and which is violated by the data under consideration ($\sum_\alpha w_\alpha^* C_\alpha = 1$). In order to prove this fact, we first observe that $\vec{w}^* \cdot \vec{C}^{\text{sep}} \le 1$ for all separable data $\vec{C}^{\text{sep}}$ (indeed $\vec{w}^* \cdot \vec{0} = 0$ for the separable data $\vec{0}$, while $\vec{w}^* \cdot \vec{C} = 1$ for the nonseparable data $\vec{C}$; we conclude using the convexity of the set of separable data). Therefore, {$\vec{w}^*, (w_i^{\text{Pauli}})^*$} represent a dual feasible [Eq. (A6)] for any separable data $\vec{C}^{\text{sep}}$. As a consequence of duality, this provides a lower bound to the primal optimum, which is $(\lambda^*)_{\text{sep}} = 0$ for separable data. Hence, $\sum_\alpha w_\alpha^* C_\alpha^{\text{sep}} + \sum_i (w_i^{\text{Pauli}})^* \le 0$, namely Eq. (A7).

### 3. General formulation of the hierarchy

**Building the correlation matrix.** The core of our approach is to build a suitable correlation matrix $\Gamma$ for

the local classical variables {$n_i$} with $n_i = (x_i, y_i, z_i)$, akin to Eq. (7). The following conditions are necessarily fulfilled for all separable data (that is, data compatible with a separable state): (1) semidefinite positivity of $\Gamma$; (2) compatibility of $\Gamma$ with the available data {$C_\alpha$}; and (3) compatibility of the local variables {$n_i$} with local qubit states (i.e., $x_i^2 + y_i^2 + z_i^2 = 1$). Crucially, verifying conditions (1)–(3) can be cast in the form of a semidefinite program. Generically, the $\Gamma$ matrix is built as

$$ \Gamma = \langle \mathbf{v}^T \mathbf{v} \rangle, \qquad (A8a) $$

$$ \mathbf{v} = [f_1(\{\mathbf{n}_i\}), \dots, f_M(\{\mathbf{n}_i\})], \qquad (A8b) $$

where the $f_i$ are arbitrary polynomial functions of the local variables $n_i$. In practice, the choice of these functions influences both the tightness and the computational cost of the approach, and considerable flexibility is offered to find the best trade-off between them—with the constraint that the data {$C_\alpha$} can be expressed as linear combinations of the entries of the $\Gamma$ matrix. Throughout this work we considered the following simple choice dictated by the nature of the data we aimed at reproducing:

$$ \mathbf{v}^{(1)} = \{1\} \cup \{n_a^{(i)}; 1 \le i \le N; a \in [x, y, z]\}. \qquad (A9) $$

As we discuss below, this defines the first level of a (convergent) hierarchy.

**Expression of the semidefinite program.** For any choice of polynomial functions $f_i$ in Eq. (A8), one may parallel the procedure described in details for qubits in Sec. A2, leading a semidefinite program of the form:

$$ \min_{X \ge 0} \langle M, X \rangle \quad \text{such that,} \qquad (A10a) $$

$$ (\text{data}) \qquad \langle A_{\alpha}^{\text{data}}, X \rangle = C_{\alpha} \qquad (\alpha \in [n_{\text{data}}]), \qquad (A10b) $$

$$ (\text{Pauli}) \qquad \langle A_i^{\text{Pauli}}, X \rangle = b_i \qquad (i \in [n_{\text{Pauli}}]), \qquad (A10c) $$

where the matrices $X, M, A$'s are as in Eq. (A3). The expression of the matrices $\Sigma_\alpha^{\text{data}}, \Sigma_i^{\text{Pauli}}$, and of the parameters $b_i$, depend on the specific choice of the monomials $v_i$. Condition (A10b) ensures that certain linear combinations of entries of the correlation matrix $\Gamma$ reproduce the (noisy) data $(1 - \lambda)C_\alpha$ (it is possible that several independent linear combinations reproduce the same data; hence, the number $n_{\text{data}}$ of such constraints may be larger than the total number $R$ of data $C_\alpha$). Conditions (A10c) enforce the compatibility of the local variables {$n_i$} with local qubit states.

**Dual problem and entanglement witness.** Exactly as for the case discussed in Sec. A2, if the (primal) problem, Eq. (A10), is unfeasible, an optimal $\lambda^* > 0$ is obtained. One may then derive an entanglement witness by solving
the (dual) problem:

$$
\begin{align}
\max_{\vec{w}, \vec{w}^{\text{Pauli}}} & \quad \vec{w} \cdot \vec{C} + \sum_{i} b_i w_i^{\text{Pauli}} \quad \text{such that,} \nonumber \\
& \quad \vec{w} \cdot \vec{C} \le 1 \nonumber \\
& \quad \sum_{\alpha} w_{\alpha} \sum_{\alpha}^{\text{data}} + \sum_{i} w_i^{\text{Pauli}} \sum_{i}^{\text{Pauli}} \le 0. \tag{A11}
\end{align}
$$

If strong duality holds (see Sec. A2), the optimal coef-
ficients $\vec{w}^*$ allow one to build an entanglement witness
as in Eq. (A7). The violation of this witness by the data
under consideration ultimately certifies the presence of
entanglement in the system.

*Convergence of the hierarchy.* Each choice of polyno-
mial functions $f_i$ in Eq. (A8) defines a different relaxation
to the set of separable data. One may actually formulate a
systematic hierarchy of such choices, converging towards
the exact separable set. The $l$th relaxation level is defined
as

$$
\begin{align}
\mathbf{v}^{(0)} &= \{1\}, \notag \\
\mathbf{u}^{(l)} &= \left\{ \prod_{r=1}^{l} n_{a_r}^{(ir)} ; a_r \in \{x, y, z\}; 1 \le i_1 \le \dots \le i_l \le N \right\}, \notag \\
\mathbf{v}^{(l)} &= \mathbf{v}^{(l-1)} \cup \mathbf{u}^{(l)}. \tag{A12}
\end{align}
$$

Throughout this work, we considered only the first relax-
ation level defined by Eq. (A9). Crucially, the correspond-
ing hierarchy converges, in the limit $l \to \infty$, 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 con-
vergence of the hierarchy [38]. It follows that the set of
correlations that can be recovered as moments of an overall
distribution $p(\{\mathbf{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 sep-
arable state, they will be detected as entangled at a finite
level of the hierarchy—although the computation cost of
high-level tests quickly increases with $l$, as one needs
to manipulate a correlation matrix of size approximately
$(4N)^l$ in the semidefinite program.

Lastly, notice that one can straightforwardly define some
hybrid levels of the hierarchy, where the entries of the vec-
tor **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 Sec. II of the main text, we already noticed that the

level-1 relaxation is left invariant by the partial transposi-
tion 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ᵢ variables [cf. Eq. (7)],
or in terms of the choice of monomials at a given relax-
ation 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 rep-
resents a hierarchy of conditions, which are completely
independent of the PPT-based criteria [33,34,39,40].

APPENDIX B: PROOF OF THE BIPARTITE ENTANGLEMENT WITNESS

Here we prove the validity of Eq. (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, \dots, N-2\}$ and $B = \{1, 3, 5, \dots, N-1\}$. We introduce local phases $\phi_a(i)$ for $a \in \{X, Y, Z\}$ and $i \in [N]$, and define

$$
2W_a = \sum_{j \in A} \sum_{j' \in B} K_{j-j'} C_{jj'}^{aa} e^{i[\phi_a(j) - \phi_a(j')]} + \text{c.c.} \quad (B1)
$$

The coefficients $K_r$ are given by

$$
K_r = K_{-r} = \frac{2}{N} \sum_{k=-\frac{N}{4}+1}^{\frac{N}{4}-1} \exp\left(\frac{2i\pi}{N}kr\right) \quad (B2)
$$

$$
= \frac{2}{N} \left[ \frac{\sin(\pi r/2)}{\tan(\pi r / N)} - \cos\left(\frac{\pi r}{2}\right) \right]. \quad (B3)
$$

The violated witnesses are then of the form

$$
W_X + W_Y + W_Z \geq -\frac{N}{2}. \qquad (B4)
$$

*Proof.* The proof of this inequality is as follows. Assuming
that the state is fully separable, we have

$$
\begin{equation}
\begin{split}
2W_a = \sum_{j,j'=0}^{N-1} & K_{j-j'} \langle a_j a_{j'} \rangle e^{i[\phi_a(j)-\phi_a(j')]} \\
& - \sum_{j,j' \in A} K_{j-j'} \langle a_j a_{j'} \rangle e^{i[\phi_a(j)-\phi_a(j')]} \\
& - \sum_{j,j' \in B} K_{j-j'} \langle a_j a_{j'} \rangle e^{i[\phi_a(j)-\phi_a(j')]}.
\end{split}
\tag{B5}
\end{equation}
$$
We use then

$$ \sum_{j,j'=0}^{N-1} K_{j-j'} a_j a_{j'} e^{i[\phi_a(j)-\phi_a(j')]} = \frac{2}{N} \sum_{k=-\frac{N}{4}+1}^{\frac{N}{4}-1} \left| \sum_{r=0}^{N-1} a_r e^{2i\pi kr/N + i\phi_a(r)} \right|^2 \ge 0. \quad (B6) $$

We then observe that for $j, j' \in 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 \neq j'$. Therefore, we have

$$ \sum_{j, j' \in A} K_{j-j'} a_j a_{j'} e^{i[\phi_a(j)-\phi_a(j')]} = \sum_{r \in A} a_r^2 - \frac{2}{N} \left| \sum_{r \in A} (-1)^{r/2} a_r e^{i\phi_a(r)} \right|^2 \le \sum_{r \in A} a_r^2. \quad (B7) $$

By the same argument, observing the $j-j'$ is an even integer for $j, j' \in B$, we have

$$ \sum_{j, j' \in B} K_{j-j'} a_j a_{j'} e^{i[\phi_a(j)-\phi_a(j')]} \le \sum_{r \in B} a_r^2. \quad (B8) $$

We therefore have

$$ 2W_a \ge -\sum_{r=0}^{N-1} \langle a_r^2 \rangle. \quad (B9) $$

Combining these inequalities for $a \in \{X, Y, Z\}$, we conclude that

$$ 2(W_X + W_Y + W_Z) \ge -\sum_{r=0}^{N-1} \langle x_r^2 + y_r^2 + z_r^2 \rangle \ge -N, \quad (B10) $$

where in the last step we use 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.

### 1. Recovering the covariance matrix criterion

In this subsection, we show how the so-called covariance matrix criterion (CMC) [26] 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_i^X, C_i^Y, C_i^Z)$, and $C_{ij}$ the $3 \times 3$ matrix with entries $C_{ij}^{ab}$ [see Eq. (1) for the definition of $C_i^a$ and $C_{ij}^{ab}$]. The CMC states that if the three-qubit state is fully separable, then there exist three real symmetric $3 \times 3$ matrices $\rho_i \ge 0$, with $\text{Tr}(\rho_i) = 1$, such that

$$ \begin{pmatrix} \rho_1 & C_{12} & C_{13} \\ C_{21} & \rho_2 & C_{23} \\ C_{31} & C_{32} & \rho_3 \end{pmatrix} \ge \begin{pmatrix} C_1^T \\ C_2^T \\ C_3^T \end{pmatrix} (C_1 \ C_2 \ C_3). \quad (C1) $$

On the other hand, the first level of our hierarchy [see Eq. (7)] implies that

$$ \Gamma = \begin{pmatrix} 1 & C_1 & C_2 & C_3 \\ C_1^T & \sigma_1 & C_{12} & C_{13} \\ C_2^T & C_{21} & \sigma_2 & C_{23} \\ C_3^T & C_{31} & C_{32} & \sigma_3 \end{pmatrix} \ge 0, \quad (C2) $$

with $\sigma_i = \begin{pmatrix} \langle x_i^2 \rangle & \langle x_i y_i \rangle & \langle x_i z_i \rangle \\ \langle x_i y_i \rangle & \langle y_i^2 \rangle & \langle y_i z_i \rangle \\ \langle x_i z_i \rangle & \langle y_i z_i \rangle & 1 - \langle x_i^2 \rangle - \langle y_i^2 \rangle \end{pmatrix}$, where $\langle \dots \rangle$ denotes an average over the (classical) $p$ distribution defining a separable state (Sec. II). Clearly, $\sigma_i \ge 0$ is a symmetric matrix of unit trace. Furthermore, positivity of $\Gamma$ 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.

### 2. Recovering the generalized spin-squeezing inequalities

Here, we show how the generalized spin-squeezing inequalities derived in Ref. [27] 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 second moments averaged over all permutations: $\{m_a, C_{aa}; a \in [x, y, z]\}$ with $m_a = N^{-1} \sum_{i=1}^N C_i^a$ and $C_{aa} = [N(N-1)]^{-1} \sum_{i \neq j} C_{ij}^{aa}$.
They consist of the following eight inequalities, valid for all fully separable states [Eq. (50) of Ref. [27]]:

$$
\begin{align*}
C_{xx} + C_{yy} + C_{zz} &\le 1, \\
C_{xx} + C_{yy} + Nm_z^2 - (N-1)C_{zz} &\le 1, \\
C_{yy} + C_{zz} + Nm_x^2 - (N-1)C_{xx} &\le 1, \\
C_{zz} + C_{xx} + Nm_y^2 - (N-1)C_{yy} &\le 1, \\
C_{xx} + N(m_y^2 + m_z^2) - (N-1)(C_{yy} + C_{zz}) &\le 1, \\
C_{yy} + N(m_z^2 + m_y^2) - (N-1)(C_{zz} + C_{xx}) &\le 1, \\
C_{zz} + N(m_x^2 + m_y^2) - (N-1)(C_{xx} + C_{yy}) &\le 1, \\
N(m_x^2 + m_y^2 + m_z^2) - (N-1)(C_{xx} + C_{yy} + C_{zz}) &\le 1.
\end{align*}
$$

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 $\mathbf{u} = (1, x_1, \sum_{i=1}^N x_i)$, and build the corresponding correlation matrix $\Gamma_x = \langle \mathbf{u}^T \mathbf{u} \rangle$, where the average is over the $p$ distribution. Using the invariance of $p$ under permutations of the qubits, we obtain

$$ \Gamma_x = \begin{pmatrix} 1 & m_x & Nm_x \\ \cdot & \langle x_1^2 \rangle & \langle x_1^2 \rangle + (N-1)C_{xx} \\ \cdot & \cdot & N\langle x_1^2 \rangle + N(N-1)C_{xx} \end{pmatrix} \ge 0. \quad (C3) $$

Positivity of $\Gamma_x$ implies that

$$ \det \begin{pmatrix} 1 & Nm_x \\ Nm_x & N\langle x_1^2 \rangle + N(N-1)C_{xx} \end{pmatrix} \ge 0 $$

and

$$ \det \begin{pmatrix} \langle x_1^2 \rangle & \langle x_1^2 \rangle + (N-1)C_{xx} \\ \langle x_1^2 \rangle + (N-1)C_{xx} & N\langle x_1^2 \rangle + N(N-1)C_{xx} \end{pmatrix} \ge 0. $$

From these inequalities, we obtain the conditions

$$ C_{xx} \le \langle x_1^2 \rangle, \quad (C4) $$

$$ Nm_x^2 - (N-1)C_{xx} \le \langle x_1^2 \rangle. \quad (C5) $$

For the same reason, we also have

$$ C_{yy} \le \langle y_1^2 \rangle, \quad (C6) $$

$$ Nm_y^2 - (N-1)C_{yy} \le \langle y_1^2 \rangle, \quad (C7) $$

$$ C_{zz} \le \langle z_1^2 \rangle, \quad (C8) $$

$$ Nm_z^2 - (N-1)C_{zz} \le \langle z_1^2 \rangle. \quad (C9) $$

Combining these inequalities, and using the property $\langle x_1^2 + y_1^2 + z_1^2 \rangle \le 1$, we recover the eight inequalities of Ref. [27]. 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. [70].

[1] R. Horodecki, P. Horodecki, M. Horodecki, and K. Horodecki, Quantum entanglement, Rev. Mod. Phys. 81, 865 (2009).

[2] J. Preskill, Quantum computing in the NISQ era and beyond, Quantum 2, 79 (2018).

[3] I. M. Georgescu, S. Ashhab, and F. Nori, Quantum simulation, Rev. Mod. Phys. 86, 153 (2014).

[4] L. Pezzè, A. Smerzi, M. K. Oberthaler, R. Schmied, and P. Treutlein, Quantum metrology with nonclassical states of atomic ensembles, Rev. Mod. Phys. 90, 035005 (2018).

[5] The necessity of entanglement to offer a quantum advantage in quantum metrology has been questioned [78]; concerning quantum simulation and computation, many efficient classical algorithms to simulate restricted families of entangled states have been developed, hence, entanglement is certainly not a sufficient ingredient to offer a quantum advantage. For pure states, however, it is necessary, as slightly entangled pure-state computations have an efficient classical representation [79].

[6] T. Monz, P. Schindler, J. T. Barreiro, M. Chwalla, D. Nigg, W. A. Coish, M. Harlander, W. Hänsel, M. Hennrich, and R. Blatt, 14-Qubit Entanglement: Creation and Coherence, Phys. Rev. Lett. 106, 130506 (2011).

[7] R. Islam, R. Ma, P. M. Preiss, M. Eric Tai, A. Lukin, M. Rispoli, and M. Greiner, Measuring entanglement entropy in a quantum many-body system, Nature 528, 77 (2015).

[8] A. M. Kaufman, M. E. Tai, A. Lukin, M. Rispoli, R. Schittko, P. M. Preiss, and M. Greiner, Quantum thermalization through entanglement in an isolated many-body system, Science 353, 794 (2016).

[9] C. Song, K. Xu, W. Liu, C.-p. Yang, S.-B. Zheng, H. Deng, Q. Xie, K. Huang, Q. Guo, L. Zhang, P. Zhang, D. Xu, D. Zheng, X. Zhu, H. Wang, Y.-A. Chen, C.-Y. Lu, S. Han, and J.-W. Pan, 10-Qubit Entanglement and Parallel Logic Operations with a Superconducting Circuit, Phys. Rev. Lett. 119, 180511 (2017).

[10] X.-L. Wang, Y.-H. Luo, H.-L. Huang, M.-C. Chen, Z.-E. Su, C. Liu, C. Chen, W. Li, Y.-Q. Fang, X. Jiang, J. Zhang, L. Li, N.-L. Liu, C.-Y. Lu, and J.-W. Pan, 18-Qubit Entanglement with Six Photons' Three Degrees of Freedom, Phys. Rev. Lett. 120, 260502 (2018).

[11] N. Friis, O. Marty, C. Maier, C. Hempel, M. Holzäpfel, P. Jurcevic, M. B. Plenio, M. Huber, C. Roos, and R. Blatt, et al., Observation of Entangled States of a Fully Controlled 20-Qubit System, Phys. Rev. X 8, 021012 (2018).
[12] A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, J. Cui, M. Rossignolo, P. Rembold, S. Montangero, T. Calarco, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Generation and manipulation of Schrödinger cat states in Rydberg atom arrays, Science 365, 570 (2019).

[13] T. Brydges, A. Elben, P. Jurcevic, B. Vermersch, C. Maier, B. P. Lanyon, P. Zoller, R. Blatt, and C. F. Roos, Probing Rényi entanglement entropy via randomized measurements, Science 364, 260 (2019).

[14] K. J. Satzinger, et al., Realizing topologically ordered states on a quantum processor, ArXiv:2104.01180 2021.

[15] F. Arute, et al., Quantum supremacy using a programmable superconducting processor, Nature 574, 505 (2019). number: 7779 Publisher: Nature Publishing Group

[16] C. S. Chiu, G. Ji, A. Bohrdt, M. Xu, M. Knap, E. Demler, F. Grusdt, M. Greiner, and D. Greif, String patterns in the doped Hubbard model, Science 365, 251 (2019).

[17] J. Koepsell, J. Vijayan, P. Sompet, F. Grusdt, T. A. Hilker, E. Demler, G. Salomon, I. Bloch, and C. Gross, Imaging magnetic polarons in the doped fermi-Hubbard model, Nature 572, 358 (2019).

[18] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator, ArXiv:2012.12281,2020.

[19] H.-S. Zhong, et al., Quantum computational advantage using photons, Science 370, 1460 (2020).

[20] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Programmable quantum simulation of 2d antiferromagnets with hundreds of rydberg atoms, ArXiv:2012.12268, 2020.

[21] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, and W. W. Ho, et al., Controlling quantum many-body dynamics in driven rydberg atom arrays, Science 371, 1355 (2021).

[22] M. Paris and J. Rehacek, *Quantum State Estimation*, Lecture Notes in Physics (Springer Berlin Heidelberg, 2004).

[23] H.-Y. Huang, R. Kueng, and J. Preskill, Predicting many properties of a quantum system from very few measurements, Nat. Phys. 16, 1050 (2020).

[24] F. Bohnet-Waldraff, D. Braun, and O. Giraud, Entanglement and the truncated moment problem, Phys. Rev. A 96, 032312 (2017).

[25] G. Blekherman, J. Gouveia, J. Helton, I. Klep, S. McCullough, J. Nie, P. Parrilo, M. Putinar, P. Rostalski, B. Sturm-fels, and R. Thomas, SIAM, edited by G. Blekherman, P. A. Parrilo, and R. R. Thomas, MOS-SIAM Series on Optimization, Vol. 13 (OUP Oxford, 2012).

[26] O. Gittsovich, P. Hyllus, and O. Gühne, Multiparticle covariance matrices and the impossibility of detecting graph-state entanglement with two-particle correlations, Phys. Rev. A 82, 032306 (2010).

[27] G. Tóth, C. Knapp, O. Gühne, and H. J. Briegel, Spin squeezing and entanglement, Phys. Rev. A 79, 042334 (2009).

[28] J. I. de Vicente and M. Huber, Multipartite entanglement detection from correlation tensors, Phys. Rev. A 84, 062306 (2011).

[29] M. Li, J. Wang, S.-M. Fei, and X. Li-Jost, Quantum separability criteria for arbitrary-dimensional multipartite states, Phys. Rev. A 89, 022325 (2014).

[30] G. Sarbicki, G. Scala, and D. Chruściński, Family of multipartite separability criteria based on a correlation tensor, Phys. Rev. A 101, 012341 (2020).

[31] J. Sperling and W. Vogel, Multipartite Entanglement Witnesses, Phys. Rev. Lett. 111, 110503 (2013).

[32] S. Gerke, W. Vogel, and J. Sperling, Numerical Construction of Multipartite Entanglement Witnesses, Phys. Rev. X 8, 031047 (2018).

[33] A. Neven, J. Carrasco, V. Vitale, C. Kokail, A. Elben, M. Dalmonte, P. Calabrese, P. Zoller, B. Vermersch, R. Kueng, and B. Kraus, Symmetry-resolved entanglement detection using partial transpose moments, ArXiv:2103.074432021.

[34] X.-D. Yu, S. Imai, and O. Gühne, Optimal entanglement certification from moments of the partial transpose, ArXiv:2103.06897,2021.

[35] M. Navascues, F. Baccari, and A. Acín, Entanglement marginal problems, ArXiv:2006.09064,2020.

[36] I. Frérot and T. Roscilde, Optimal entanglement witnesses: a scalable data-driven approach, ArXiv:2101.02038,2021.

[37] F. Baccari, D. Cavalcanti, P. Wittek, and A. Acín, Efficient Device-Independent Entanglement Detection for Multipartite Systems, Phys. Rev. X 7, 021042 (2017).

[38] J. B. Lasserre, Global optimization with polynomials and the problem of Moments, SIAM J. Optim. 11, 796 (2001).

[39] A. Peres, Separability Criterion for Density Matrices, Phys. Rev. Lett. 77, 1413 (1996).

[40] A. C. Doherty, P. A. Parrilo, and F. M. Spedalieri, Detecting multipartite entanglement, Phys. Rev. A 71, 032333 (2005).

[41] R. F. Werner, Quantum states with einstein-podolsky-Rosen correlations admitting a hidden-variable model, Phys. Rev. A 40, 4277 (1989).

[42] W. K. Wootters, Entanglement of Formation of an Arbitrary State of Two Qubits, Phys. Rev. Lett. 80, 2245 (1998).

[43] C. Gross and I. Bloch, Quantum simulations with ultracold atoms in optical lattices, Science 357, 995 (2017).

[44] T. Fukuhara, S. Hild, J. Zeiher, P. Schauß, I. Bloch, M. Endres, and C. Gross, Spatially Resolved Detection of a Spin-Entanglement Wave in a Bose-Hubbard Chain, Phys. Rev. Lett. 115, 035302 (2015).

[45] H. Levine, A. Keesling, A. Omran, H. Bernien, S. Schwartz, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, High-Fidelity Control and Entanglement of Rydberg-Atom Qubits, Phys. Rev. Lett. 121, 123603 (2018).

[46] V. Lienhard, S. de Léséleuc, D. Barredo, T. Lahaye, A. Browaeys, M. Schuler, L.-P. Henry, and A. M. Läuchli, Observing the Space- and Time-Dependent Growth of Correlations in Dynamically Tuned Synthetic Ising Models with Antiferromagnetic Interactions, Phys. Rev. X 8, 021070 (2018).

[47] A. Keesling, A. Omran, H. Levine, H. Bernien, H. Pichler, S. Choi, R. Samajdar, S. Schwartz, P. Silvi, and S. Sachdev, et al., Quantum kibble-Zurek mechanism and
critical dynamics on a programmable rydberg simulator, Nature 568, 207 (2019).

[48] S. de Léséleuc, V. Lienhard, P. Scholl, D. Barredo, S. Weber, N. Lang, H. P. Büchler, T. Lahaye, and A. Browaeys, Observation of a symmetry-protected topological phase of interacting bosons with rydberg atoms, Science 365, 775 (2019).

[49] A. Browaeys and T. Lahaye, Many-body physics with individually controlled rydberg atoms, Nat. Phys. 16, 132 (2020). number: 2 Publisher: Nature Publishing Group

[50] P. Fersterer, A. Safavi-Naini, B. Zhu, L. Gabardos, S. Lepoutre, L. Vernac, B. Laburthe-Tolra, P. B. Blakie, and A. M. Rey, Dynamics of an itinerant spin-3 atomic dipolar gas in an optical lattice, Phys. Rev. A 100, 033609 (2019).

[51] A. Patscheider, B. Zhu, L. Chomaz, D. Petter, S. Baier, A.-M. Rey, F. Ferlaino, and M. J. Mark, Controlling dipolar exchange interactions in a dense three-dimensional array of large-spin fermions, Phys. Rev. Res. 2, 023050 (2020).

[52] L. Gabardos, B. Zhu, S. Lepoutre, A. M. Rey, B. Laburthe-Tolra, and L. Vernac, Relaxation of the Collective Magnetization of a Dense 3D Array of Interacting Dipolar S=3 Atoms, Phys. Rev. Lett. 125, 143401 (2020).

[53] M. W. Doherty, N. B. Manson, P. Delaney, F. Jelezko, J. Wrachtrup, and L. C. L. Hollenberg, The nitrogen-vacancy colour centre in diamond, Phys. Rep. 528, 1 (2013).

[54] J. Choi, H. Zhou, S. Choi, R. Landig, W. W. Ho, J. Isoya, F. Jelezko, S. Onoda, H. Sumiya, and D. A. Abanin, et al., Probing Quantum Thermalization of a Disordered Dipolar Spin Ensemble with Discrete Time-Crystalline Order, Phys. Rev. Lett. 122, 00 (2019).

[55] P. Jurcevic, B. P. Lanyon, P. Hauke, C. Hempel, P. Zoller, R. Blatt, and C. F. Roos, Quasiparticle engineering and entanglement propagation in a quantum many-body system, Nature 511, 202 (2014).

[56] P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith, M. Foss-Feig, S. Michalakis, A. V. Gorshkov, and C. Monroe, Non-local propagation of correlations in quantum systems with long-range interactions, Nature 511, 198 (2014).

[57] J. Zhang, G. Pagano, P. W. Hess, A. Kyprianidis, P. Becker, H. Kaplan, A. V. Gorshkov, Z.-X. Gong, and C. Monroe, Observation of a many-body dynamical phase transition with a 53-qubit quantum simulator, Nature 551, 601 (2017). number: 7682 Publisher: Nature Publishing Group

[58] K. A. Landsman, C. Figgatt, T. Schuster, N. M. Linke, B. Yoshida, N. Y. Yao, and C. Monroe, Verified quantum information scrambling, Nature 567, 61 (2019).

[59] M. K. Joshi, A. Elben, B. Vermersch, T. Brydges, C. Maier, P. Zoller, R. Blatt, and C. F. Roos, Quantum Information Scrambling in a Trapped-Ion Quantum Simulator with Tunable Range Interactions, Phys. Rev. Lett. 124, 240505 (2020).

[60] C. Monroe, W. Campbell, L.-M. Duan, Z.-X. Gong, A. Gorshkov, P. Hess, R. Islam, K. Kim, N. Linke, and G. Pagano, et al., Programmable quantum simulations of spin systems with trapped ions, Rev. Mod. Phys. 93, 00 (2021).

[61] P. Wittek, Algorithm 950, ACM Trans. Math. Softw. 41, 21:1 (2015).

[62] Available at http://www.mosek.com/.

[63] https://github.com/ifrerot/SDP-multipartite-entanglement.

[64] L. Mazza, D. Rossini, R. Fazio, and M. Endres, Detecting two-site spin-entanglement in many-body systems with local particle-number fluctuations, New J. Phys. 17, 013015 (2015).

[65] I. Frérot, P. Naldesi, and T. Roscilde, Entanglement and fluctuations in the XXZ model with power-law interactions, Phys. Rev. B 95, 245111 (2017).

[66] I. Frérot, P. Naldesi, and T. Roscilde, Multispeed Prethermalization in Quantum Spin Models with Power-Law Decaying Interactions, Phys. Rev. Lett. 120, 050401 (2018).

[67] P. Krammer, H. Kampermann, D. Bruß, R. A. Bertlmann, L. C. Kwek, and C. Macchiello, Multipartite Entanglement Detection via Structure Factors, Phys. Rev. Lett. 103, 100502 (2009).

[68] M. Cramer, M. B. Plenio, and H. Wunderlich, Measuring Entanglement in Condensed Matter Systems, Phys. Rev. Lett. 106, 020401 (2011).

[69] P. Hauke, R. J. Sewell, M. W. Mitchell, and M. Lewenstein, Quantum control of spin correlations in ultracold lattice gases, Phys. Rev. A 87, 021601 (2013).

[70] G. Vitagliano, P. Hyllus, I. n. L. Egusquiza, and G. Tóth, Spin Squeezing Inequalities for Arbitrary Spin, Phys. Rev. Lett. 107, 240502 (2011).

[71] S. Sachdev, *Quantum Phase Transitions* (Cambridge University Press, Cambridge, UK, 2011).

[72] I. Frérot and T. Roscilde, Quantum Critical Metrology, Phys. Rev. Lett. 121, 020402 (2018).

[73] P. Hauke, M. Heyl, L. Tagliacozzo, and P. Zoller, Measuring multipartite entanglement through dynamic susceptibilities, Nat. Phys. 12, 778 (2016).

[74] G. Mathew, S. L. L. Silva, A. Jain, A. Mohan, D. T. Adroja, V. G. Sakai, C. V. Tomy, A. Banerjee, R. Goreti, A. V. N. R. Singh, and D. Jaiswal-Nagar, Experimental realization of multipartite entanglement via quantum fisher information in a uniform antiferromagnetic quantum spin chain, Phys. Rev. Res. 2, 043329 (2020).

[75] A. Scheie, P. Laurell, A. M. Samarakoon, B. Lake, S. E. Nagler, G. E. Granroth, S. Okamoto, G. Alvarez, and D. A. Tennant, Witnessing entanglement in quantum magnets using neutron scattering, ArXiv:2102.08376, 2021.

[76] O. Gühne and G. Tóth, Entanglement detection, Phys. Rep. 474, 1 (2009).

[77] X.-G. Wen, Choreographed entanglement dances: Topological states of quantum matter, Science 363, eaal3099 (2019).

[78] D. Braun, G. Adesso, F. Benatti, R. Floreanini, U. Marzolino, M. W. Mitchell, and S. Pirandola, Quantum-enhanced measurements without entanglement, Rev. Mod. Phys. 90, 035006 (2018).

[79] G. Vidal, Efficient Classical Simulation of Slightly Entangled Quantum Computations, Phys. Rev. Lett. 91, 147902 (2003).