Inferring non-linear many-body Bell’s inequalities from average two-body correlations: Systematic approach for arbitrary spin-j ensembles

,


I. INTRODUCTION
Multipartite entanglement is a central feature of quantum many-body systems, fundamentally challenging our ability to efficiently simulate them on classical computers [1,2].For the same reason, quantum entanglement distributed among many degrees of freedom represents a key resource for quantum simulators and computers.Consequently, proving that the multipartite states prepared in quantum simulators or computers are indeed entangled -namely, the task of entanglement certification -is a key step in assessing the quantum advantage offered by such devices.Depending on the assumptions made about the individual components of the device, two different paradigms may appear suitable.In a socalled device-dependent framework, the subsystems are well characterized: the Hilbert space is known (e.g. a qubit space), and the measurements correspond to welldefined quantum observables (e.g.spin measurements); in this framework, entanglement certification relies on the violation of a certain entanglement witness [3].On the other hand, in a device-independent framework, no assumption is made about the Hilbert space of the subsystems, and consequently the measurements correspond to unknown quantum observables; this framework appears especially suitable when considering effective fewlevel systems, where the actual Hilbert space can contain an unlimited number of physical degrees of freedom.* irenee.frerot@gmail.com Relaxing certain assumptions about the system clearly makes entanglement certification more demanding; nevertheless, device-independent entanglement certification is possible if the violation of a certain Bell's inequality [4] can be established.Designing many-body Bell tests is the focus of the present paper.
As fully characterizing the many-body correlations among the subsystems requires exponentially-many measurements, any scalable Bell test must rely on incomplete information, obtained from an accessible number of measurements -for instance, the knowledge of few-body correlations among the parties [5][6][7].In particular, in order to mitigate scalability issues, a successful strategy is to symmetrize the data from which Bell non-locality is to be certified, over all permutations of the parties [7][8][9][10][11][12].Correspondingly, the Bell's inequalities relevant to such coarse-grain features of the system involve a number of coefficients which is independent of the number of parties.On the experimental side [9,10], this allows one to reduce the statistical uncertainty on the data entering the Bell's inequality; on the theoretical side, this allows one to reduce the computational complexity, often leading to the analytical characterization of Bell's inequalities [7,8,12].A second challenge for entanglement certification is to take advantage, to the largest possible extent, of all the available information, without a priori knowing the structure of entanglement within the many-body state.In particular, failing to violate all known Bell's inequalities does not imply that device-independent entanglement certification is impossible based on the available data.This motivates the development of data-driven methods [5,7,11], where the data serve as input into FIG.1: Illustration of an ideal multipartite Bell test for entanglement certification.Right: A composite quantum system prepared by the source is shared among N spatially-separated observers (on the sketch, N = 3), each of which chooses a measurement settings a ∈ {0, 1, . . ., k − 1}, obtaining an outcome s ∈ {−j, −j + 1, . . ., j}.In this work, we focus on spin measurements on quantum spin-j particles, sketched as Stern-Gerlach magnets oriented along directions na -but Bell tests are independent of these assumptions.This procedure is repeated several times in order to accumulate statistics, where at each round the measurement settings and the observed outcomes may vary.If the observed statistics exhibit Bell's non-local correlations, one certifies that the multipartite state is quantum-mechanically entangled.Left: In a Bell test, each subsystem i is treated as a black box, with no assumption about the Hilbert space.The actual observables ŝ(i) a being measured are attributed arbitrary input labels a, and the measurement outcomes s have no physical meaning.Entanglement is therefore certified in a device-independent manner.
an algorithm which builds, from the data themselves, a tailored Bell's inequality.In the case of two-outcome measurements -especially suited to spin-1/2 systems -a data-driven method for permutationally-invariant Bell's inequalities based on semi-definite-positive relaxations, has been proposed in the past [11].To our knowledge, this method has however not lead to the discovery of new families of Bell's inequalities; and its extension to more outcomes -relevant to spin-j > 1/2 systems -has never been achieved.
Taking inspiration from statistical physics, in this work we develop an alternative flexible data-driven method which takes, as input data, one-and two-body correlation functions averaged over all permutations of the subsystems -for any number of measurement outcomes and settings.Similarly to the method of ref. [11], the complexity of our algorithm is independent of the system size, and tests exhaustively an infinite number of Bell's inequalities in a data-driven fashion.This lead us to recover tigher versions of all previously known permutationallyinvariant Bell's inequalities [7,8,12] in an unbiased way.Furthermore, in contrast to ref. [11], our scheme is directly applicable to any number of outcomes, and is validated by the study of quantum spin j > 1/2 ensembles.Finally, the Bell's inequalities inferred by our method are non-linear in the input data, and tightly wrap around the polytope of local-variable models (see Fig. 2 for an example).This feature offers a fundamental scaling improvement regarding experimental requirements, including for all previously-known Bell's inequalities invariant under permutations [7][8][9][10]12].Among other results obtained with our novel method, we discover new families of many-body Bell's inequalities, for measurements involving arbitrarily-many outcomes, violated by paradigmatic many-body entangled states for ensembles of quantum spins j ≥ 1/2 -namely, spin singlets and spin-squeezed states -a topic of timely experimental revelance to many experimental platforms manipulating qudit ensembles.
Before entering into the details of our new method, in the rest of this section we review the framework of deviceindependent entanglement certification (I A), and the notion of Bell's inequalities invariant under permutations (I B).In Section II, we present our method in the case of two-outcomes measurements (II A), and apply it to improve over and extend previously-known Bell's inequalities in the case of spin-singlets (II C) and spin-squeezed states (II D) for spin-1/2 ensembles.In particular, we emphasize the fundamental scaling improvement offered by the non-linear nature of the Bell's inequalities inferred by our algorithm (see e.g.Fig. 2).Section III is then devoted to the hitherto-unexplored case of spin-j > 1/2 (namely, qudits) ensembles, for which we extend our approach to an arbitrary number of outcomes (III A).We then apply it again to spin-singlets (III B) and spin-squeezed state (III C), leading us to characterize analytically novel families of Bell's inequalities, valid for any number of parties and outcomes.Section IV contains experimental considerations: in Section IV A, we list different platforms and their respective capabilities to detect entanglement and Bell correlations; in Section IV B we discuss the statistical requirements to acquire the data used as input to our algorithm.Finally, Section V displays our conclusions and prospects.
A. Device-independent entanglement certification Bell scenario.Arguably, the violation of Bell's inequalities [4] represents the most robust scheme to certify entanglement, avoiding detailed assumptions about the physical nature of the degrees of freedom being measured, and about the accurate calibration of the measurements being performed.In this so-called device-independent scenario (Fig. 1), each subsystem i ∈ {1, . . .N } (for instance, a quantum spin-j) is treated as a black box, namely, no assumption is made about the actual Hilbert space of the system.This black-box treatment is especially relevant when dealing with effective few-level systems.On each subsystem, k different measurement settings can be implemented.In practice, they correspond to certain quantum operators ŝ(i) a (a ∈ {0, . . ., k −1}), for instance spin measurements along particular directions n a , but in a device-independent scenario the actual quantum measurement which is performed is not assumed; instead, only the outcome of the measurement, denoted s, is collected (see Fig. 1) (throughout the paper, we denote as Ô a quantum observable, and as O the outcome of its measurement).The only assumption made is that the number d of possible outcomes for each measurement is finite.In practice, the possible values of s are the eigenvalues of the quantum operator ŝ(i) a (for instance, the 2j + 1 possible values of a spin-j measurement), but in a device-independent scenario these are mere labels for the outcomes, with no specific physical meaning.For convenience and later connection with quantum violations of Bell's inequalities when performing appropriate spin measurements on quantum many-body systems, we denote these possible outcomes as s ∈ {−j, −j + 1, . . ., j} with d = 2j +1 -but it should be emphasized that within a device-independent framework, these labels are arbitrary.A Bell experiment [4] consists in repeating the following sequence: 1) choose a setting a (i) ∈ {0, . . ., k − 1} for each subsystem; 2) perform the corresponding measurements, yielding the N outcomes s = {s (i) } N i=1 .By repeating this sequence, varying the measurement settings a = {a (i) } N i=1 , statistics of the measurement outcomes are collected.Complete information is obtained if one reconstructs all N -body marginal probability distributions: P (s|a) for all choices of settings a.If one denotes Π(i) a,s the projector onto the eigenspace of the observable ŝ(i) a associated to the eigenvalue s, then these probabilities are given by where ρ is the density-matrix of the system -notice that even if, in a device-independent scenario, we remain agnostic about the Hilbert space over which ρ acts, such a decomposition exists in principle.
Bell's inequalities and entanglement certification.The state ρ is not entangled (namely, it is separable) if it can be decomposed as a statistical mixture of product states: where ρ(i) λ is an arbitrary local quantum state (pure or mixed) for subsystem i, acting on the local Hilbert space whose dimension is arbitrary.λ is some classical random variable, sampled with probability measure dµ(λ), which encodes classical correlations among the local quantum states ρ(i) λ .The central observation behind device-independent entanglement certification is that if ρ is separable [Eq.( 2)], then P (s|a) can always be reproduced by a local-variable (LV) model in the sense of J. S. Bell [4,13]: where In a device-independent framework, we do not know the explicit expressions of the projectors Π(i) a,s corresponding to the measurements which are actually being performed; and even the Hilbert space of the system, over which the quantum state ρ and these projectors are defined, is unknown and arbitrary -it could even be infinite dimensional.Yet, regardless of the actual Hilbert space describing the system, if the state is not entangled, then a decomposition as in Eq. ( 3) must exist for the experimentally-observed correlations contained in P (s|a).
Therefore, if conversely P (s|a) is found to violate a Bell's inequality -denying the possibility to decompose P (s|a) as in Eq. (3)-, then ρ must be entangled.Crucially, this holds regardless of the Hilbert space of the individual subsystems, and regardless of the measurements which were actually performed to generate P (s|a).Violating a Bell's inequality therefore certifies that ρ is entangled in a device-independent manner.
Note that in principle, violating a Bell's inequality allows for quantum information protocols more powerful than merely witnessing entanglement [4], which is the task on which we focus in this paper.

B. Permutationally-invariant Bell's inequalities from two-body correlations
Certifying entanglement from two-body correlations.
Overall, reconstructing P (s|a) requires collecting k N (d N − 1) probabilities.This exponential scaling clearly makes full reconstruction of these marginals unpractical, and therefore, methods based on partial information have been developed.The simplest non-trivial strategy, which we follow in this paper, is to consider jointly all two-body marginals: P (ij) (s, t|a, b) (i = j), namely the probability to obtain the pair of outcomes (s, t) if measurement a is performed on subsystem i, and measurement b on subsystem j, for all possible pairs of subsystems 1 ≤ i < j ≤ N , and all possible pairs of measurement settings 0 ≤ a, b ≤ k−1.
Local-variable models as distributions over classical spin configurations.The à-la-Bell formulation of LV models as in Eq. ( 3) makes transparent the link with entangled quantum states.It is however more intuitive to represent LV models as probability distributions over the measurement results treated as classical variables s (i) a [7].Indeed, as first proved by Fine [4,14], a LV decomposition as in Eq. ( 3) exists if and only if there exists a grand-probability distribution P LV [σ] over the fictitious ensemble of classical variables σ = {s (i) a ; a = 0, . . ., k − 1; i = 1, . . ., N }, such that the observed statistical properties are obtained as marginals against P LV , i.e. [7,14]: where δ is the Kronecker symbol (δ x,y = 1 if x = y, and δ x,y = 0 otherwise).In LV models, measurement results may therefore be viewed as sampled from an underlying classical "spin" configuration σ, where k "hidden" d−level spins s (i) a ∈ {−j, −j + 1, . . .j} are attached to each subsystem i, encoding the outcome of the measurement.While, at each measurement run with setting a = {a (i) } N i=1 , the value of only one of the k hidden spins is revealed [namely, s (i) ] also objectively exist independently of the act of their measurement.This contradicts standard interpretations of quantum physics if they correspond to incompatible quantum observables performed on the same subsystem, [ŝ b ] = 0 -and is categorically excluded if the P (ij) (s, t|a, b) violate a Bell's inequality, and if actionsat-a-distance are not allowed [4].
Permutationally-invariant Bell's inequalities.Deciding whether or not the marginals P (ij) (s, t|a, b) are compatible with a grand-probability P LV (σ) can be mapped onto a so-called inverse Ising problem [7], which can generically be solved in polynomial time by Monte-Carlo methods -while worst-case instances are exponentially hard.A convergent hierarchy of relaxations to this problem has also been developed, whose computational cost is strictly polynomial at each relaxation level [5].Here, we drastically simplify the problem by further symmetrizing the data over all permutations of the subsystems [8,11], which leads us to introduce: Bell's inequalities are constraints, of the form: where f is some function, and B c the so-called classical bound, obeyed by all distributions PLV (s, t|a, b) which descend from a grand-probability P LV (σ).If the particular P under investigation happens to violate such a Bell's inequality [namely, if f ( P ) < B c ], then no grandprobability P LV (σ) can ever explain the data, which in turn implies that the quantum state ρ of the system must be entangled.
Our main result is to construct a very flexible datadriven algorithm, whose complexity is independent of N , allowing one to build a Bell's inequality violated by the data P (s, t|a, b) (Section II A).This allows us to recover all previously known permutationally-invariant Bell's inequalities which are robustly violated by appropriate quantum states in the thermodynamic limit [7,9], to improve these Bell's inequalities by considering more measurement settings (Section II), and to generalize them to scenarios with arbitrarily-many outcomes (Section III).We discuss the potentialities of several experimental platforms to observe Bell non-locality in Section IV, and draw our conclusions in Section V.

II. TWO-OUTCOMES MEASUREMENTS
Summary of the main results.In this section, we introduce our method by focusing on the simplest situation where the measurements can only deliver d = 2 outcomes.The method itself is presented in Section II A: the key results are contained in Eqs. ( 9), ( 10) and (11), which form the core of our data-driven algorithm.To be practically useful, the algorithm must be fed with carefully-chosen quantum data.In Section II B, we consider a situation where the data correspond to spin measurements on a collection of quantum spin-1/2.We expose the dependence of these data on collective-spin fluctuations, as represented by Eq. (14).We then begin our data-driven exploration of Bell's inequalities with spin singlets (Section II C) and spin-squeezed states (Section II D), for which permutationally-invariant Bell's inequalities are already known.In both cases, we find tighter Bell's inequalities, leading to sufficient "witness" conditions on collective spin fluctuations which are easier to satisfy than existing ones.Concerning singlets (Section II C), our main finding is a family of Bell's inequalities for arbitrarily-many measurement settings [Eq.(15)].The corresponding witness condition is contained in Eq. (17).Concerning spin-squeezed states (Section II D), we illustrate the generic improvement offered by the non-linear nature of our Bell's inequalities on Fig. 2. We then go beyond existing Bell's inequalities [8][9][10] by adding an extra measurement setting, whose advantage for entanglement certification is illustrated on Figs. 3 and 4.

A. A convex-optimization algorithm
We assume that all measurements ŝ(i) a can only deliver d = 2 possible outcomes, denoted s = ±1/2 (the usual convention in the literature would be to denote them s = ±1, but we follow our general convention s ∈ {−j, −j + 1, . . .j} with d = 2j + 1; as already emphasized, these labels are arbitrary).Instead of working with the pair probability distribution P (ij) (s, t|a, b), we equivalently consider one-and two-body correlations s (the two representations are related by elementary linear transformations).As coarse-grain features of the experimental data, equivalently to the averaged pair distribution P (s, t|a, b), we consider the one-and twobody correlations summed over all permutations of the subsystems: In a LV description, the s a are N k classical Ising spins (with values ±1/2).A LV model compatible with the (coarse-grain) experimental data corresponds to a probability distribution P LV ({s (i) a }) over the configurations of these Ising spins, such that M a and C ab are obtained as marginals against P LV .Let us assume that a LV model fitting the data exists, and derive necessary conditions obeyed by the corresponding M a and C ab (namely, Bell's inequalities).We first introduce the collective variables a , and their fluctuations δS a = S a − S a , so that we have: The terms on the r.h.s of Eq. (8b) are not directly observable.In particular, the terms s correspond to correlations among the measurement settings a and b on the same subsystem i.In the general case, these settings correspond to incompatible quantum observables, [ŝ b ] = 0, and therefore these terms do not have a direct meaning in quantum physics.However, they are perfectly well-defined in LV models.The first key observation, which underlies the method developed in the present paper, is that for any k × k positive semi-definite (PSD) matrix A 0, and for any configuration of the collective variables S a , we have a,b δS a A ab δS b = δS T AδS ≥ 0. We introduced the vector notation S := (S 0 , . . .S k−1 ) T , and used the fact that, by definition of a PSD matrix, u T Au ≥ 0 for any vector u.Therefore, for any A 0 and any vector h = (h 0 , . . .h k−1 ) T , we have: where . This is a Bell's inequality, obeyed by all data (C ab , M a ) compatible with LV models, any PSD matrix A, and any vector h.The bound E max (A, h) may easily be evaluated by enumerating all 2 k configurations of the s variables, whenever k (the number of settings) is not too large.The goal is then to find a PSD matrix A and a vector h such that Eq. ( 9) is violated.In order to build them, our second key observation is that E max (A, h) is a convex function of its arguments.A simple proof of convexity, inspired by statistical physics, is to write and to recognize that log Z β is a convex function for any β.Furthermore, Tr(A C) + h • M, which is a linear function of A and h, is also convex.Therefore, we may introduce the convex cost function: which by Eq. ( 9) is non-negative if (C, M) are compatible with a LV model.Our data-driven algorithm [15] consists therefore in solving the following optimization problem: As the PSD constraint A 0 maintains the convex nature of the optimization problem [16], if there exists a Bell's inequality of the form of Eq. ( 9) which is violated by the data, then we have the guarantee to find the corresponding A 0 and h s.t.L(A, h) < 0. Notice that if L(A, h) = −l < 0 in Eq. ( 10) then for any x > 0, L(xA, xh) = −xl, so that L is unbounded below.In a practical implementation of the algorithm, one may therefore add a cutoff on A and h; for this work, we have imposed Clearly, the possibility to discover new and useful Bell's inequalities via our method crucially depends on the input quantum data {C ab , M a }, which must be able to display Bell's non-locality.We will consider a situation where the quantum data are obtained by spin measurements (Sec.II B).In this case, C ab and M a are completely determined by the first-and second-moments of the collective spin.As a first application, we will recover and improve over the existing Bell's inequalities in scenarios with d = 2 outcomes, which are violated by appropriate measurements on spin singlets [7] (Sec.II C) and spinsqueezed states [9] (Sec.II D), and whose violation is robust in the thermodynamic limit.In Section III, we will then generalize these results to scenarios with d > 2 outcomes.

B. Spin measurements
Throughout the paper, we investigate the violation of Bell's inequalities when the local measurement settings correspond to spin measurements in the xy-plane, in a direction independent of the subsystem.We emphasize that this choice is only a convenient way to produce hypothetical quantum data, used as input to our data-driven algorithm.The discovered Bell's inequalities themselves are valid independently of any assumption about the system.Furthermore, we present in details several Bell's inequalities discovered by our algorithm; this lead us to derive simple conditions on the quantum state of a spin ensemble which are sufficient to violate the Bell's inequalities if the appropriate measurements are performed (in the literature, such conditions are often referred to as "Bell-correlation witnesses").These witness conditions are independent of the specific data we used to discover the Bell's inequalities of interest.We choose therefore: where Ŝ(i) x and Ŝ(i) y are local spin observables in directions x and y.
ŝ(i) a defines a projective spin measurement along the direction (cos θ a , sin θ a ), and has therefore eigenvalues ±1/2.Introducing the collective spin Ĵx = y , we also define the collective spin observables: With these conventions, the quantum data used as input of our algorithm, and against which the Bell's inequalities are evaluated, are (for a given quantum state ρ): where δ Ĵa = Ĵa − M a , so that δ Ĵa δ Ĵb = Ĵa Ĵb + Ĵb Ĵa /2 − M a M b is the covariance of the collective spin observables Ĵa and Ĵb .We defined Â = Tr( Âρ), and we introduce the variance Var( Â) = Â2 − Â 2 .

C. A family of Bell's inequalities for singlet-like correlations
As a first application of our data-driven method, we derive Bell's inequalities maximally violated by many-body singlets, defined by Ĵ2 x = Ĵ2 y = Ĵ2 x = 0. Many-body singlets are zero-eigenstates of the total spin operator Ĵ2 = Ĵ2 x + Ĵ2 y + Ĵ2 z .They are therefore SU (2) invariant (that is, they are left invariant by any rotation exp[−in • Ĵ] with n a unit vector) and generalize the Bell [17]-all entangled -, and are naturally produced as ground states of Heisenberg antiferromagnets [18], e.g. at low energy in Fermi-Hubbard models [19].We emphasize that the working assumption of having a many-body singlet is only used to produce ideal quantum data, which then serve as input to our algorithm, leading us to discover new Bell's inequalities.The Bell's inequalities themselves, and the corresponding Bell correlation witnesses, are independent of any assumption about the quantum state.It is already known that a state is entangled when Var( Ĵx ) + Var( Ĵy ) < N/4 [20].It is also known that [ Ĵ2 x + Ĵ2 y ]/N < 1/(8 + 6 √ 2) ≈ 0.060660 . .., which is a more demanding condition, leads to violation of a manybody Bell's inequality [7].The measurement strategy to maximally violate the Bell's inequality of ref. [7] is composed of k coplanar spin measurements at angles θ a = aπ/k.Our main result in this Section is to show that the Bell's inequality of ref. [7] is not the tightest one in this measurement scenario for k ≥ 4, leading us to discover a new family of Bell's inequalities.We find that Bell-nonlocality can be demonstrated whenever [Var( Ĵx ) + Var( Ĵy )]/N < 1/2 − 4/π 2 ≈ 0.094715, in the limit of k → ∞.
Bell's inequality.As input quantum data, we consider a perfect spin singlet, for which M a = 0 and 4 14), and using the property Ĵa ρsinglet = 0 for any collective spin operator Ĵa and any singlet state ρsinglet ].Applying our algorithm to these data for up to k = 10, we find that the following Bell's inequality is violated: where on the second line we used Eq. ( 9).The classical bound B c is obtained by noting that Quantum violation.To evaluate Eq. ( 15) against a generic quantum state [Eq.( 14)], not necessarily SU (2) invariant, we first introduce the matrix (16) Violation of the Bell's inequality is therefore detected whenever: The tightest condition is achieved in the limit k → ∞, yielding the bound 1/2−4/π 2 .This condition requires no assumption about the underlying quantum state -apart from being composed of N individual spin-1/2-, and only assumes a correct calibration of the measurements of Ĵx and Ĵy .Notice that the condition of Eq. ( 17) is tighter than the one derived in ref. [7], not only because the r.h.s is larger -and therefore detecting more data as exhibiting Bell's non-locality in the same measurement scenario -, but also because it involves variances of the collective operators, making the condition more robust against experimental noise -see the related Fig. 2, and Section IV B.
State-of-the-art Bell's inequality.In the context of Bell's non-locality, spin-squeezing is known to be essential for the robust violation of the following Bell's inequality involving k = 2 measurement settings per subsystem [8][9][10]22]: Notice that we are using the convention that the outcomes are ±1/2, different from the convention ±1 used in the above-cited papers; this explains why the coefficients of Eq. ( 18) are different.We have used the experimental ref.[9] to infer data serving as input to our algorithm, leading us to recover a tighter version of the above Bell's inequality: This shows that, if Eq. ( 18) had not been known from ref. [8], the tighter Eq. ( 19) would have been recovered in a data-driven way by our method.This clearly demonstrates the concrete advantage offered by our method in analyzing experimental data in an unbiased way.Notice that while the coefficients of the Bell's inequalities Eq. ( 18) and Eq. ( 19) are the same, Eq. ( 18) involves the correlations C ab [Eq.(7b)], and not Cab = C ab − M a M b [Eq.(8b)].Therefore, Eq. ( 19) includes the extra term −(M 0 − M 1 ) 2 ≤ 0, and is therefore strictly tighter than Eq. ( 18) -see Fig. 2 for an illustration.Since this extra term is of order O(N 2 ) while the classical bound is O(N ), the relative improvement generically grows with N .The classical bound is found, following Eq.( 9), by writing B = (δS , and noting that (s 0 − s 1 ) 2 + s 0 + s 1 ≥ −1 for all possibles values of s 0 , s 1 = ±1/2.This Bell's inequality can be violated by preparing a spinsqueezed state, defined by N Var( Ĵy ) < Ĵx 2 [21], and performing two projective spin measurements in directions ŝ(i) x cos θ ± Ŝ(i) y sin θ [9,22].Computing the quantum value from Eq. ( 14), we obtain B = 4 sin 2 θVar( Ĵy ) − 2 cos θ J x − N sin 2 θ.The optimal angle θ, minimizing B for fixed data (Var( Ĵy ), Ĵx ), is cos θ = Ĵx /[N − 4Var( Ĵy )].For this choice of measurements, we obtain Notice that for Eq. ( 18), a similar condition may be derived, but involving Ĵ2 y instead of Var( Ĵy ).Whenever Ĵy = N with = 0, Ĵ2 y ∼ N 2 , which represents a fundamental obstruction to the violation of Eq. ( 18) in the thermodynamic limit for non-ideal data.Instead, working with the tighter Eq. ( 19), and the corresponding criterion involving Var( Ĵy ), such obstruction is removed.For perfect squeezed states [Var( Ĵy ) → 0, Ĵx → N/2], we can obtain violation up to B = −5N/4.In this Section, we show that the robustness of Bell non-locality detection for spin-squeezed states can be improved by considering extra measurements (k ≥ 3).
Finding tightest and more robust Bell's inequalities.To find better Bell's inequalities, our strategy was to consider quantum data [Eq.( 14)] obtained from a squeezed state at the limit of violating the Bell's inequality Eq. ( 19), and add extra measurements in the xy-plane to potentially discover other violated Bell's inequalities.In particular, adding a third spin measurement Ŝ(i) y along the y-axis, we found a family of Bell's inequalities, defined by the following coefficients [see Eq. ( 9)]: where a ≥ 0. The corresponding Bell's inequality reads: and reduces to Eq. ( 19) when a = 0. Remarkably, for a = 1, Eq. ( 21) represents a tighter version of a Bell's  19), which involves k = 2 measurement settings at angles ±θ.
Dashed-dotted-dotted green line: relative violation of the Bell's inequality of Eq. ( 21), which involves a third measurement along the y-axis (see the sketch on the bottomleft corner), with a = 1 in Eq. ( 21).Solid orange line: same inequality, for the optimal value of the parameter a.
The quantum state, chosen from ref. [9], has a mean spin 2 Ĵx /N = mx = 0.98, and transverse collective-spin fluctuations 4Var( Ĵy)/N = χ 2 = 0.272 (assuming Ĵy = 0).For all inequalities, the absolute value of the relative violation is equal to the amount of white noise tolerated by the data to observe a non-zero violation (see text).
inequality analyzed in ref. [12] -tighter, due to the nonlinear nature of Eq. ( 21) which involves C instead of C. Similarly to Eq. ( 19), we discovered this family of Bell's inequalities parametrized by a in a data-agnostic way, using data from a spin-squeezed state as input to our algorithm.The classical bound B c may be found in the following manner.Noting that x T Ax = (x 0 − x 1 + ax 2 ) 2 for any vector x, we may write: The classical bound B c = −N (1 + a/2) 2 is then found by enumerating the configurations of the variables s (i) a = ±1/2.On the other hand, in the quantum measurement setting we consider, we find: The quantum data ( Cab , M a ) being fixed, it is then natural to consider the optimal values of θ and a to have the most robust violation of the Bell's inequality.If white noise is added to the data, then (C ab , M a ) → (1 − r)(C ab , M a ) with r the noise amplitude.If we assume that Ĵy = 0, so that Var( Ĵy ) = Ĵ2 y , then correspondingly B → (1 − r) B .The noise robustness may be interpreted as the intrinsic robustness of a given Bell's inequality violation against generic errors during the preparation of the quantum system, modelled as ρ = (1 − r)ρ (ideal) + r1/D with D the dimension of the total Hilbert space (a more detailed discussion of the experimental requirements to accurately estimate the data is given in Section IV B).Maximizing the noise robustness is therefore equivalent to maximizing the ratio , where we introduced m x = 2 Ĵx /N (the Rabi contrast [9]) and χ 2 = 4Var( Ĵy )/N the scaled second moment.For each value of θ, we may then find the value of a which maximizes this ratio.As illustrated on Fig. 3 for the data of ref. [9] (m x = 0.98 and χ 2 = 0.272 assuming that Ĵy = 0), adding a third measurement along the y-axis and optimizing over the parameter a yields a systematic improvement over both Eq. ( 19) (which involves only two measurement settings), and over Eq. ( 21) with a = 1, as proposed in ref. [12].Notice that we assumed that Ĵy = 0. Therefore, in the specific measurement settings we considered, working with the C quantities in the Bell's inequalities Eq. ( 19) or Eq. ( 21) is equivalent to C; however, in analyzing experimental data where Ĵy is never exactly zero, the non-linear nature of our Bell's inequalities (namely, working with C) will lead to a systematic improvement over all previously known Bell's inequalities [8][9][10]12].It would be interesting to analyze the experimental data of refs.[9,10] from this perspective -this would certainly lead to a more robust detection of Bell correlations.
Finally, for given values of (m x , χ 2 ), we may find the measurement angle θ and parameter a which maximize the robustness of the violation of Eq. ( 21).This is done on Fig. 4, which shows the parameter regime in the (m x , χ 2 ) where Bell non-locality is detected.For comparison we also plot the regime where non-locality is detected based on the violation of Eq. ( 19), on the violation of Eq. ( 21) with a = 1, and where entanglement is detected based on the Wineland spin squeezing criterion m 2 x > χ 2 .The Bell's inequality Eq. ( 21) with optimal a systematically extends the parameter space where non-locality can be detected with k = 3 settings.Notice that this parameter space can be further extended by considering more measurement settings [12].
Further improvement.Even more robust Bell's inequalities may be found by considering k ≥ 2 pairs of measurements in the xy plane, in directions ŝ(i) x cos θ a − Ŝ(i) y sin θ a , and one measurement along y denoted ŝ(i) k = Ŝ(i) y .In this Bell scenario with 2k+1 settings applied to spin-squeezed states, this generically leads to Bell' inequalities of the  21) with the optimal parameter a; dasheddotted-dotted green line: same Bell's inequality with a = 1 [12]; dashed red line: violation of the Bell's inequality Eq. ( 19) [9].Black star: experimental data from ref. [9], assuming Ĵy = 0. Notice that the previously known results were in fact involving Ĵ2 y .We have proved in this work that they remain valid with Var( Ĵy) instead of Ĵ2 y , leading to systematically tighter criteria.In particular, the data of ref. [9] (black star), using the actual value of Ĵy in the experiment, would be lower along the vertical axis, and similarly for the data of ref. [10] (not shown).(24) for some data-tailored coefficients α a and β a .Similarly as in the previous paragraph, for given values of (m x , χ 2 ), it is then possible to numerically optimize over the measurement angles θ a in order to maximize the violation robustness.Similarly to the case with three measurements discussed above, one can expect to obtain systematically tighter Bell's inequalities as compared to ref. [12], where an arbitrary number of settings are considered in a similar measurement scenario.

III. ARBITRARY-OUTCOMES MEASUREMENTS
Summary of the main results.In the previous section, we focused on measurements with d = 2 outcomes -and considered a physical implementation with spin-1/2 measurements.In this section, we extend these results to arbitrarily-many outcomes (d > 2), corresponding to the physical situation where spin measurements are performed on individual spin-j components (with d = 2j + 1).Sec. III A presents an incremental generalization of the algorithm of Sec.II A, which incorporates an extra feature of the quantum data [Eq.(25)].This turns out to be an essential ingredient to generalize the Bell's inequalities of the previous section.We first consider spin-j singlets in Sec.III B, and restrict our attention to k = 3 spin measurements in a given plane.We unveil an increasingly complex situation for j > 1/2, as illustrated on Fig. 5, with 2j inequivalent Bell's inequalities, already in this simple setting.We could however characterize analytically two families of Bell's inequalities which emerged from our algorithm (Sec.III B).One of them extends Eq. ( 15) to arbitrary half-integer spins (for k = 3 measurements), and the corresponding witness condition is given by Eq. (39).The other family is valid for both integer and half-integer spins, and the witness condition is given by Eq. ( 41).We conclude our exploration in Sec.II D with spin-j squeezed states.Our main result is a generalization of the Bell's inequality of Eq. ( 19) to arbitrary j ≥ 1/2 [Eq.( 43)].The corresponding witness condition for spin-j squeezed states is given by Eq. (45).i =j P (ij) (s, t|a, b) may then be obtained by averaging over all permutations of the subsystems.In Appendix A, we give a general formulation of our datadriven algorithm for finding a Bell's inequality violated by P .However, aiming at finding new Bell's inequalities for many-spin systems with j > 1/2, we found sufficient to include only: as an extra coarse-grain feature of the quantum data, in addition to M a and C ab defined in Eq. (7).Apart from this modification, we may then follow the same construction as in Section II A, where the kN local classical variables s a can now take the d possible values −j, −j + 1, . . .j.The analogue of Eq. ( 9) now contains an extra term h (2) • M (2) to allow for Bell's inequalities involving this extra feature of the data.Explicitly, for any PSD matrix A, and any vectors h = (h 1 , . . .h k ) and h (2) = (h k ), we have: where now E max (A, h, h (2) ) = max s∈{−j,...j} k E(s), with a s 2 a ].Eq. ( 26) is a Bell's inequality, satisfied by all data M a , M (2) a and C ab compatible with a LV model with d-outcome measurements.We may then parallel the end of Section II A: introduce the convex cost function L(A, h, h (2) ) = Tr(A C)+h•M+h (2) •M (2) +N E max (A, h, h (2) ), and minimize it via a convex-optimization routine, imposing the PSD constraint A 0. If we find L < 0, a violated Bell's inequalities is then reconstructed from the corresponding A, h and h (2) .
Applying this algorithm, we discovered Bell's inequalities violated by spin-j spin singlets, and by spin-j squeezed states.These Bell's inequalities generalize the results of Section II to arbitrary j ≥ 1/2.

B. Bell's inequalities for arbitrary-j many-body singlets
We start our investigation of Bell's inequalities tailored to spin-j systems by considering, as input to our algorithm, many-body singlets.This will lead us to extend some the results obtained in Section II C for j = 1/2 to arbitrary spins.Here again, the assumption of having a many-body singlet is only used to produce quantum data leading us to discover new Bell's inequalities via our data-driven algorithm.The Bell's inequalities are independent of any assumption about the quantum systems being measured, and the witness inequalities only assume that a collection of N spin-j particles are measured along appropriately-calibrated axes.As in the case of j = 1/2, spin singlets are SU (2)-invariant states defined by the sole condition Ĵ2 x = Ĵ2 y = Ĵ2 z = 0.It is known that if Var( Ĵx ) + Var( Ĵy ) + Var( Ĵz ) ≤ N j, then the state is multipartite entangled [23], and therefore spin singlets are entangled for any j.We have considered k = 3 coplanar spin measurements, in directions ŝ(i) a = Ŝ(i) x cos θ a + Ŝ(i) y sin θ a , with {θ a } = (t 1 , 0, −t 2 ).We did not find violated Bell's inequalities with k = 2 settings, and using non-coplanar spin measurements did not lead to more robust Bell's inequalities.However, we do not exclude that better Bell's inequalities could be found using non-coplanar measurements, with k ≥ 4 settings, or including more general SU (2j + 1) measurements.In summary, this setting appeared as the simplest one to discover new Bell's inequalities, and even in this simplest scenario we could not characterize all the Bell's inequalities which appear when increasing j.Fig. 5 summarizes our findings, where we plot the violation of Bell's inequalities in the (t 1 , t 2 ) plane, for j ∈ {1/2, 1, 3/2, 2}, together with the witness condition on the collective spin variance Var( Ĵx ) to observe violation (assuming global SU (2) invariance).We shall discuss in details two families of Bell's inequalities, which were characterized analytically for arbitrary j.

General considerations
We begin with general considerations on the Bell's inequalities discovered by using, as input to our algorithm, the data obtained measuring a many-body singlet along k directions in the xy-plane.As a consequence of SU (2) invariance, a many-body singlets has no mean spin orientation: for any direction a, Ĵa = a ] 2 = N j(j + 1)/3.In a many-body singlet, the quantum data Cab , M a and M (2) a are [Eqs.( 8) and ( 25)]: (in the specific cases discussed in further details below, θ 0 = t 1 , θ 1 = 0 and θ 2 = −t 2 ).
Structure of the Bell's inequalities.As a consequence of M a = 0, the Bell's inequalities tailored to singlets do not involve terms linear in M a [namely, in the notations of Eq. ( 26), we have h = 0], and they take the general form: where A is a symmetric PSD matrix.For LV models, we have: where we defined Ãab = A ab − h a δ ab , and: Notice that the bound is tight for N even.Indeed, if {s a=0 is also saturating E max .We may therefore  (27)] to our algorithm, we considered k = 3 spin measurements in the xy-plane, forming angles t1, 0, and −t2 with the x axis (inset of panel a).We found violated Bell's inequalities as in Eq. ( 28), with a conventional a ] 2 = 1, and a classical bound given by Eq. ( 30) . (a,b,c,d): Difference between the classical bound and the quantum value [Eq.( 31)], divided by N .(a) j = 1/2; (b) j = 1; (c) j = 3/2; (d) j = 2.The integers 1−5 label inequivalent families of Bell's inequalities, which are found in the different regions of parameters (t1, t2) (boundaries between these regions are indicatively emphasized as dotted white lines).In general, we find 2j Bell's inequalities inequivalent under relabelling of the outcomes.Labels 1 and 2 respectively correspond to the coefficients of Eqs.(38) and (40).(e) Along t := t1 = t2 [white dashed line on panels (ad)], bound on Var( Ĵx)/N to observe Bell non-locality assuming SU (2) invariance in this measurement setting [Eq.(37)].The maximum at t = π/3 for half-integer spins corresponds to label 1, and the witness condition is that of Eq. ( 39).The right-most local maximum at t = arccos[1/(4j)] corresponds to label 2, and the witness condition is that of Eq. (41).Quantum value on a singlet.Considering the quantum value on a spin singlet [Eq.( 27)], we find: Using that Ãab = Ãba , we have ab Ãab cos(θ a − θ b ) = ab e −iθa Ãab e iθ b .The optimal angles, leading to the maximal violation of the Bell's inequality, are those which maximize ab e −iθ b Ãab e iθa .We therefore define: This should be contrasted to the case of LV models where, in order to find the classical bound, one maximizes ab s a Ãab s b over the variables s a ∈ {−j, . . .j} [see Eq. ( 30)].On Fig. 5(a-d), we plot the quantum violation of the Bell's inequalities we found, for j ∈ {1/2, 1, 3/2, 2}, with k = 3 spin-measurement directions {θ a } = (t 1 , 0, −t 2 ).In general, varying the measurement angles t 1 and t 2 , we find 2j inequivalent inequalities.We characterized analytically one Bell's inequality appearing for all half-integer j, and one family appearing for all j (see below).
Witness condition.For each Bell's inequality of the form Eq. ( 28) found by our approach, and for given measurement directions {θ a }, one may derive witness conditions which can be measured via global measurements on an ensemble of N spin-j particles, and which demonstrate the capability of the quantum state to violate the considered Bell's inequality without further assumptions (in particular, without assuming SU (2) invariance).We first express the average value of the Bell operator [Eqs.(28) and (29)] in terms of spin observables: where δ Ĵa = Ĵa − Ĵa , Ĵa = N i=1 ŝ(i) a , and ŝ(i) y sin θ a .Using A ab = A ba , Ãab = Ãba and elementary algebra, we obtain: where {x, ŷ} = xŷ + ŷx, and: We derive explicitly the corresponding witnesses for two families of Bell's inequalities in the next Section.For the sake of illustration, on Fig. 5(e) we plot the witness condition for the Bell's inequalities we found for j ∈ {1/2, 1, 3/2, 2}, with angles {θ a } = (−t, 0, t) [that is, along the diagonal of panels (a-d) of the same Fig.5].
To realize this plot and derive a simple condition involving only the collective spin fluctuations Var( Ĵx ), we have further assumed SU (2) invariance of the state.With this assumption, we find: where B singlet is given by Eq. ( 31).The witness condition is B SU (2) < −N E max ( Ã) with E max ( Ã) given by Eq. ( 30).Explicitly, the witness condition for SU (2)invariant states is: where the min is over {s a } ∈ {−j, . . ., j} k .This upper bound is plotted on Fig. 5

(e).
2. A Bell's inequality for half-integer spin singlets.The Bell's inequality presented in Eq. ( 15), and valid for j = 1/2, can be extended to arbitrary spins.Here, we focus on the simplest extension with k = 3 measurement settings, which corresponds to the region labelled '1' on panels (a) and (c) of Fig. 5.The coefficients we found are: We notice that for arbitrary complex numbers x 0 , x 1 , x 2 , we have: For integer spins, choosing s 0 = s 1 = s 2 = 0 gives E max = 0, so that the classical bound cannot be violated by measuring spin singlets [Eq.(31)].In contrast, for half-integer spins, the minimal value of |s 0 − s 1 + s 2 | is 1/2, so that E max = −1/4.Instead, for singlets, choosing the measurement angles {θ a } = (π/3, 0, −π/3) (which is optimal), we find ab e −iθa Ãab e iθ b = 0. We will now derive a witness condition for this optimal choice of measurements, valid with no assumption about the quantum state.From Eq. (35), we find that a , we find B = (9/2)[Var( Ĵx ) + Var( Ĵy )].Therefore, violation of the Bell's inequality Eq. ( 28), with coefficients given by Eq. ( 38), and whose classical bound is B c = N/4, is possible whenever: In order to reach this conclusion, we only assumed that a collection of N spin-j particless, with j a half-integer, is measured along well-calibrated axes x and y.Maximal violation is obtained for perfect singlets which satisfy Var( Ĵx ) = Var( Ĵy ) = 0.This generalizes a result of ref. [7] to arbitrary half-integer spins.On Fig. 5(e) where SU (2) invariance is further assumed (so that Var( Ĵx ) = Var( Ĵy )), this condition corresponds to the maximum at t = π/3 for j ∈ {1/2, 3/2}.
3. A family of Bell's inequalities for arbitrary spin singlets.
Further improvement.The Bell's inequalities reported here are the simplest ones discovered via our datadriven method, involving only k = 3 co-planar spin measurements.Adding extra measurements, possibly genuine SU (d) measurements, can only lead to more robust Bell's inequalities, and looser witness conditions than Eqs.( 39) and (41).We leave this exploration open to future studies, for which our algorithm [15] represents a natural starting point.
C. Spin-squeezed states A Bell's inequality for arbitrary-spin squeezed states with two settings.Similarly to the measurement scenario leading to the violation of Eq. ( 19) for spin-1/2 squeezed states, we consider a situation where two projective spin measurements, ŝ(i) y sin θ, are locally performed on a collection of N spinj particles.Using data from a spin-j squeezed state, we find the following generalization of Eq. ( 19): For j = 1/2, we have M = N/4, and therefore recover Eq. (19).More generally, the classical bound ≥ 0 is obtained by writing , where E(s 0 , s 1 ) = (s 0 − s 1 ) 2 − 2s 2 0 − 2s 2 1 + s 0 + s 1 = −(s 0 + s 1 )(s 0 + s 1 − 1).Since s 0/1 are the outcomes of spin-j measurements, they are either both integers, either both half-integers.This implies that s 0 + s 1 is always an integer.Since E(x) = −x(x − 1) ≤ 0 for all integers x, we conclude that B ≥ 0 for LV models.On the other hand, for the measurement setting we consider, we have the quantum value: x ] 2 .The optimal measurement angle θ, leading to the minimal value of B , is s.t.cos θ = Ĵx /[4N s 2 x − 4Var( Ĵy )] (if this is ≤ 1), for which we obtain: Violation is detected whenever B < 0. This generalizes the results for j = 1/2 [9,22] to arbitrary spins.In Appendix B, we present another Bell's inequality, for which violation has been found with j ≤ 1.Clearly, adding extra measurements (k ≥ 3) could only lead to more robust Bell's inequalities (see Sec. II D for the case j = 1/2): we leave this exploration open to future works.

IV. EXPERIMENTAL IMPLEMENTATION
In this work, we have introduced a new methodology to learn, from experimental data themselves, the best device-independent entanglement criterion that the data allow one to construct -in the form of a Bell's inequality whose coefficients are inferred via a data-driven algorithm.We have demonstrated the effectiveness of this new approach by using, as input to our algorithm, either actual experimental data [9], or data which could be obtained by collecting the appropriate two-body correlations on realistic many-body quantum states.The new Bell's inequalities presented in this work include and surpass all robust permutationally-invariant Bell's inequalities reported so far in the literature.Therefore, these Bell's inequalities could already be useful for entanglement certification in existing or near-term quantum devices, if the relevant states are prepared, and if the appropriate measurements are performed (see below).But most importantly, the data-driven nature of our approach makes it especially suitable to explore a virtually-infinite variety of experimental data, potentially unveiling new and unexpected Bell's inequalities.It is indeed not unrealistic to anticipate that present-day quantum simulators and computers are processing quantum-entangled states, while the experimentalists are not able to prove it simply because they lack entanglement criteria tailored to their experimental data.To facilitate this exploration by other researchers, we have released a pedagogical openaccess version of the code used throughout this paper [15].Furthermore, in this section, we present an (incomplete) list of experimental platforms able to produce suitable data, allowing one to potentially certify the preparation of many-body entangled state with the data-driven method presented in this work.The present section does not contain new results; rather, it consists of a guide towards the existing literature relevant to the experimental implementation of device-independent entanglement certification.Such implementation requires answering two questions: • A) Can one collect the data sets used as inputs to our data-driven algorithm?
• B) Can one prepare quantum many-body states manifesting Bell's non-locality?
The measurement question can take two conceptuallydifferent forms: A1) the one-and two-body correlations forming the data set [Eqs.(7), ( 8) and ( 25)] are measured individually, which requires individual addressing of the subsystems; A2) these data are inferred from the fluctuations of collective observables.In the first case (A1), one realizes a situation conceptually close to a proper Bell test, even though avoiding e.g. the locality loophole might be very challenging.In contrast, in the second case (A2), one does not realize a Bell test, but rather demonstrates the ability to prepare many-body entangled states which would yield violation of the reconstructed Bell's inequalities, if such a Bell test could be implemented.Clearly, (A2) is less demanding and requires only access to collective variables (as realized in [9,10]), which are sums of spins or pseudo-spins of individual components of the system.These individual components could be atoms of spin j = 1/2, 1, 3/2, . . ., or atoms, ions and superconducting qubits realizing effective few-level systems.Firstand second-moments of the collective-spin components must be measured.For spin-j particles, witness conditions such as those we derived [Eq.( 41) and ( 45)] may involve quantities of the form Ô = , where Ŝ(i) n is a spin observable along an arbitrary direction n, and f (x) may be an arbitrary function.Such observables can be measured via collective measurements in atomic systems, e.g. after Stern-Gerlach splitting of the magnetic sublevels before imaging [24].Indeed, we may express Ô = j s=−j is the number of atoms detected with spin s after the Stern-Gerlach splitting with a magnetic field gradient along n.
Concerning the preparation (B), the goal is to generate strongly entangled many-body states, such as squeezed states or spin singlets.Notice that while spin singlets and spin squeezed states are paradigmatic examples of many-body entangled states, on which we focused in this work to demonstrate the effectiveness and flexibility of our new method, other classes of states are potentially interesting; for instance, Dicke states are clear candidates [25,26] while the potentialities of j ≥ 1 ensembles subjected to more general SU (d) measurements, are virtually unlimited.In this section, we present an incomplete list of platforms for which both A and B questions can be answered positively.We then discuss the experimental effort, in terms of the number of repetitions of the preparation-measurement procedure, in order to establish the violation of permutationallyinvariant Bell's inequalities.

A. Experimental platforms
Atomic ensembles.These are clouds of (not-necessarily cold) atoms with spin.A) The total spin components can be measured employing quantum Faraday effect, i.e. looking at the polarisation rotation of the light passing through the atomic cloud [27].This method is also frequently termed as spin polarisation spectroscopy (SPS).In principle, one has access here to the full quantum statistics of the total spin components.Using standing driving fields, one can also detect spatial Fourier components of the collective spin [28].B) Atomic ensembles are particularly suitable to achieve strong squeezing of the atomic spin.Using quantum feedback, spin singlet states have also been prepared [29][30][31].Combining feedback with the ideas of Ref. [28], practically arbitrary spin-spin correlations could be generated [32].Using a one-dimensional atom-light interface, quantum spin noise limit can be achieved [33].Spin-squeezed states and even Bell's non-locality have been achieved for macroscopic ensembles, using optical-cavity-assisted measurements [10].
Ultracold spinor Bose-Einstein condensates.A) The techniques developed for atomic ensembles can be applied to Bose-Einstein condensates [34].In ultracold trapped spinor gases, the principal non-linear mechanism leading to squeezing (among other interesting entangled states) corresponds to spin-changing collisions [35].Here, all spin components and their fluctuations can be measured.Beyond spin components, e.g. the nematic tensor for spin-1 condensates can be measured, including in a spatially-resolved way [36].B) The possibility to violate the many-body Bell's inequalities of Ref. [8] were in fact first confirmed in twin-modes squeezed states in spinor condensates [9], and these systems allow one to generate non-classical states going beyond squeezing (for a review, see [21]).Entanglement between spatially separated condensates, and even Einstein-Podolsky-Rosen steering, can be generated [37][38][39].More recently, high-spin cold-atoms ensembles, which display dipolar magnetic interactions, have been generated [40][41][42][43], which appear as ideal playgrounds to explore Bell's inequalities tailored to arbitrary-spin systems, as established in the present work.
Ultracold atoms in optical lattices.A) Ultracold atoms in optical lattices provide one of the best platforms for quantum simulations [44].All the methods mentioned above can be carried over to atoms in optical lattices.SPS has emerged as a promising technique for detecting quantum phases in lattice gases via the coherent mapping of spin-correlations onto scattered light, realizing quantum non-demolition measurements.In particular, spatially-resolved SPS that employs standing wave laser configurations [28] allows for a direct probing of magnetic structure factors and order parameters [45][46][47][48].Moreover, quantum gas microscopes, which are able to resolve individual atoms located in single lattice sites, have been developed [49][50][51].These techniques allow for a direct inspection into the spatial structure of entanglement within the system [52], and direct violation of the Bell's inequalities (A1) could be envisioned.B) These systems may lead to a very large variety of strongly-correlated many-body states.In particular, spin singlets, which are ground states of quantum antiferromagnets according to theorems by Mattis [18], are expected to emerge at low-energy in quantum simulators of the Fermi-Hubbard model [19].A review of potentially achievable correlations can be found in Ref. [53].High-spin atomic ensembles in optical lattices [54,55] clearly represent very promising systems to investigate novel classes of entangled many-body states, especially concerning Bell's inequalities involving many outcomes.
Trapped ions.A) Trapped ions represent a very versatile platform, and one of the most promising candidates for quantum computing.Small systems of ions allow for full quantum tomography of the density matrix, from which the statistics of all observables can be recovered.This includes in particular spin-spin correlations (see for instance [56]), but also e.g.third order correlations which have been used for the detection of genuine threebody entanglement in a system of few ions.Clearly, trapped ions are a platform of choice to directly probe the spatial structure of quantum correlations, and direct violation of the Bell's inequalities could be achieved (A1).B) Few-ion systems can be used for the generation of a very wide variety of entangled states on demand: recent examples include ground states of lattice gauge theory models, [57][58][59], and dynamically-generated entanglement [60], among others.
Significant others.These include, but are not limited to, atoms in nano-structures [61], Rydberg atoms [62], and a large number of condensed matter systems, ranging from circuit QED, through quantum dots, to superconducting Josephson junctions [63].All of these systems could potentially prepare and detect the entangled states suitable to violate the Bell's inequalities investigated in this work -and most importantly, generate correlation patterns from which data-driven methods such as ours could reveal novel Bell's inequalities.

B. Measurement effort
In order to implement the data-driven method presented in this paper, one needs to estimate the data M a and Cab [Eq.( 8)] (the generalization to include also terms such as M (2) a [Eq.( 25)] for measurements with d > 2 outcomes is straightforward).In order to evaluate them directly to realize a device-independent entanglement test (i.e.without inferring them via collective measurements to estimate Bell correlation witnesses), one needs to have the experimental capability to individually address each subsystem, and to choose independently the measurement setting a ∈ {0, . . ., k − 1} on each of them.The following procedure may be repeated R times.
We denote as R (i) a the number of times the setting a has been implemented on subsystem i, and R (ij) ab the number of times the pair of settings (a, b) has been implemented on the pair of subsystems (i, j), i.e.: On average, for each subsystem, each setting a is implemented R/k times; and for each pair of subsystems, each pair of settings (a, b) is implemented R/k 2 times.The data are then obtained as: The collective quantity M a is typically scaling as O(N ) with fluctuations of order O( √ N ) (this holds whenever there is a finite correlation length in the system).On the other hand, the collective quantity C ab scales as O(N 2 ) with fluctuations of order O(N ).The quantity Cab is instead scaling as O(N ), but its fluctuations, stemming from the fluctuations of C ab and M a M b which are both of order O(N ), are also of order O(N ).Therefore, the error on M a and Cab due to finite statistics scale as: R , respectively.The most demanding estimation is for the two-body correlations contained in Cab .Notice that, as a consequence of the fact that the data involve only extensive quantities, the number R of measurement runs required to reach a given relative precision of scales as R ∼ k 2 / 2 , and therefore does not scale with the system size.Notice also that if the goal is not to collect the data to be used as input of our data-driven algorithm, but instead to evaluate the violation of a given permutationally-invariant Bell's inequality, such as those presented in this work, it might be more efficient to select the measurement settings with probabilities depending on the coefficients of the Bell's inequality in question.One then needs to estimate the error on ab A ab Cab + a h a C a − B c (with B c the classical bound), which should be significantly negative for the certification to be conclusive (under the gaussian-statistics assumption, see below).
Improvement due to the non-linear nature of the Bell's inequalities.As already noticed in Section II D, the fact that our Bell's inequalities involve the (non-linear) Cab quantities leads to significantly tighter results than the Bell's inequalities involving C ab (see in particular Fig. 2), especially for N large.Indeed, C ab is typically of order O(N 2 ), while M a , Cab = O(N ), together with the classical bound which is also O(N ) [see Eq. ( 9) or Eqs. ( 15) and (19) for explicit examples].Therefore, any systematic error will lead to an error of O(N 2 ) on C ab , making more challenging in practice the detection of Bell non-locality for large N (see, however, refs.[9,10]).Instead, given that all terms in our Bell's inequalities are extensive, a systematic error will lead to O(N ) deviations, which does not represent an obstruction to scalable Bell tests.
Relaxing the gaussian-statistics assumption.If the implicit gaussian-statistics assumption leading to the scaling R ∼ k 2 / 2 is to be relaxed, more elaborate finitestatistics analysis must be carried on, typically using tail bounds on the distribution of outcomes (see e.g.refs.[64][65][66]).Similarly, if a Bell correlation witness, based on collective measurements, is to be evaluated, from which the ability of the prepared state to violate a Bell's inequality is to be assessed, special care in the data analysis must be taken if the gaussian-statistics assumption is relaxed [12], leading to an overhead in terms of measurement effort.
C. Summary of the concrete implementation of our method.
In summary, detecting entanglement via our datadriven method proceeds in four steps: 1. Define a partition of the multipartite system into N subsystems, and select several (incompatible) local quantum observables ŝ(i) a for i ∈ {1, . . .N }, whose outcome are denoted s , either by measuring individually the subsystems, or by inferring such data via collective measurements; 3. Use these data as input to our algorithm to potentially find a violated Bell's inequality [15].If no violation is found, one could modify the measurements chosen at step (1); 4. Analytically analyze the Bell's inequality inferred from the data, to understand the essential features leading to entanglement detection.

V. CONCLUSIONS
We have presented a new data-driven method to detect multipartite entanglement in quantum simulators and computers.We devised an algorithm (Section II A) which constructs a violated non-linear Bell's inequality from one-and two-body correlations averaged over all permutations of the subsystems.Our approach is applicable to any number of measurement outcomes.In order to do so, we have expressed the two-body coefficients of the Bell's inequality as a positive semidefinite matrix, whose optimization allows for a systematic exploration of all potentially violated Bell's inequalities of this form.As an illustration of the potentialities of this new approach to entanglement detection, we could improve over previously-known many-body Bell's inequalities violated by j = 1/2 spin-squeezed [8][9][10]12] and spin singlet states [7] in the thermodynamic limit (Section II).In addition, we could extend these results to similar states for arbitrary j > 1/2 individual spins, by considering Bell scenarios with arbitrarily-many outcomes (Section III) -to our knowledge, this represents the first example of such families of many-body Bell's inequalities.As our (non-linear) Bell's inequalities involve only zero-momentum fluctuations, sufficient conditions on many-body quantum states for their violation could be established, in the form of Bell-correlation witnesses -involving first moments and variances of collective observables.Such witnesses can be measured in state-of-theart cold-atoms systems with only global measurements (Section IV).Importantly, the non-linear nature of the Bell's inequalities reconstructed by our method offers a fundamental scaling improvement over the linear Bell's inequalities which have been considered so far.Due to its very flexible nature and a very small computational cost (independent of the system size, and exponential in the number of measurement outcomes), our data-driven approach opens the way to the systematic exploration of permutationally-invariant Bell's inequalilites in manyqudits systems -as a matter of fact, the Bell's inequalities presented in the paper represent only a fraction of all those discovered with our approach [15], already for the simple classes of spin-squeezed and spin-singlet states.We anticipate that exploring other many-body entangled states, for instance considering Dicke states, or going beyond spin measurements to consider genuine SU (d) measurement [36], either theoretically, or directly from experimental data, will lead to the discovery of yet many other and -by construction -useful many-body Bell's inequalities.Finally, we would like to point out that our algorithm searches for Bell's inequalities whose two-body coefficients form a positive semi-definite matrix.While all robust permutationally-invariant Bell's inequality reported in the literature satisfy this condition, it is worth investigating its limitations.

FIG. 3 :
FIG.3:Bell's inequality violation for j = 1/2 spin-squeezed states, as a function of the measurement angle θ.Dashed red line: relative violation of the Bell's inequality of Eq. (19), which involves k = 2 measurement settings at angles ±θ.Dashed-dotted-dotted green line: relative violation of the Bell's inequality of Eq. (21), which involves a third measurement along the y-axis (see the sketch on the bottomleft corner), with a = 1 in Eq. (21).Solid orange line: same inequality, for the optimal value of the parameter a.The quantum state, chosen from ref.[9], has a mean spin 2 Ĵx /N = mx = 0.98, and transverse collective-spin fluctuations 4Var( Ĵy)/N = χ 2 = 0.272 (assuming Ĵy = 0).For all inequalities, the absolute value of the relative violation is equal to the amount of white noise tolerated by the data to observe a non-zero violation (see text).

1 FIG. 5 :
FIG.5:Bell's inequalities for many-body spin singlets.As input quantum data [Eqs.(27)] to our algorithm, we considered k = 3 spin measurements in the xy-plane, forming angles t1, 0, and −t2 with the x axis (inset of panel a).We found violated Bell's inequalities as in Eq. (28), with a conventional normalization ab A 2 ab + a [h while saturating the bound, by choosing the configuration {s (opt) a } k−1 a=0 for N/2 subsystems, and {−s (opt) a } k−1 a=0 for the other N/2 subsystems.