Probing quantum entanglement from magnetic-sublevels populations: beyond spin squeezing inequalities

Spin squeezing inequalities (SSI) represent a major tool to probe quantum entanglement among a collection of few-level atoms, and are based on collective spin measurements and their ﬂuctuations. Yet, for atomic ensembles of spin-j atoms and ultracold spinor gases, many experiments can image the populations in all Zeeman sublevels s = − j, − j + 1 , . . . , j , potentially revealing ﬁner features of quantum entanglement not captured by SSI. Here we present a systematic approach which exploits Zeeman-sublevel population measurements in order to construct novel entanglement criteria, and illustrate our approach on ground states of spin-1 and spin-2 Bose-Einstein condensates. Beyond these speciﬁc examples, our approach allows one to infer, in a systematic manner, the optimal permutationally-invariant entanglement witness for any given set of collective measurements in an ensemble of d - level quantum systems.


Introduction
To prepare and detect quantum-entangled states has become a central goal for experimental manybody physics. On the one hand, it demonstrates the ability to probe regimes unthinkable within a classical framework, either at equilibrium or during the dynamics. On the other hand, it represents a crucial step towards using such systems as resources for quantum-enhanced applications, such as sensing [17,31], or the simulation of quantum many-body problems [15]. Here, our main focus are those many-body systems which are most easily probed by collective measurements, such as atomic ensembles and ultracold spinor gases. In such systems, entanglement can be famously revealed by socalled spin-squeezing inequalities (SSI), which involve measurements of first-and second-moments of collective spin operators [22,29,31,36,37,39,41,42,44]. SSI and their generalizations [16,39,41,42] have found an impressively broad range of successful applications, from detecting entanglement in cold and hot atomic ensembles [23,31] to unveiling many-body Bell non-locality [8,30,32,35], from probing quantum phase transitions [11] and quantum quenches [5, 6] to applications in quantum-enhanced metrology [31]. In particular, a handful of generalized SSI is sufficient to capture the complete palette of entanglement that can be revealed via collective-spin measurement in j = 1/2 spin ensembles [39]. In the context of j = 1 spinor Bose-Einstein condensates (BEC), by restricting the spin states to effective two-dimensional subspaces, entanglement has also been probed via those SSI [7, 18,24,25,33]. In parallel, the last decade has witnessed the development of several experiments investigating much-larger-j spinor gases [1, 3, 4, 9, 10, 13, 14, 18, 19, 24-28, 33, 40, 43, 45], mixtures [2,20], or effective qudits ensembles [34]. In these systems, the populations of all Zeeman sublevels s = −j, −j + 1, . . . , j can be measured, typically by fluorescence imaging after collective spin rotations and Stern-Gerlach splitting [9, 10, 14, 18, 24-26, 28, 33, 45]. The SSI generalized in refs. [42] to spin j > 1/2 can be potentially be applied in such systems in order to detect entanglement. Yet, the entanglement patterns incorporated by j > 1/2 spin ensembles is potentially much richer than what can be captured by SSI, with new classes of entangled states and new entanglement mechanisms. The optimal use of such collective imaging data to probe entanglement beyond SSI has remained an open problem, for which our paper offers a novel approach.
In order to do so, we actually formulate and solve a more general problem: considering a set of single-atom observables, we assume that their average value and pair correlations, averaged over all permutations of the atoms, can be estimated. In ref. [41], a similar framework was proposed, and new entanglement witnesses were derived. Yet, in order to potentially find one violated inequality, in general the latter framework requires to test a number of inequalities which grows exponentially with the number of local observables, such as populations in different Zeeman sublevels as we consider here. In this work, we instead derive a single inequality based on such data to reveal entanglement. This inequality not only summarizes all known SSI [39,41,42], but also detects new forms of entanglement. Specializing then to Zeeman-sublevels population measurements in spinor gases, we apply this general approach to discover novel families of entanglement criteria, akin to but not captured by SSI. We illustrate our findings on representative states of a spin-1 BEC, as produced and characterized in several experiments [9, 10, 18,33,45]. Finally, we extend these findings to the largely unexplored case of a spin-2 BEC.
Framework. Our goal is to exclude the possibility to decomposeρ, the many-body density operator of the system, as a statistical mixture of product states over individual atoms:ρ =ρ sep := λ p λ ⊗ N i=1ρ λ (i). Here,ρ λ (i) is an arbitrary internal state (pure or mixed) for atom i = 1, . . . , N , and p λ > 0 is an arbitrary probability distribution, defining a separable stateρ sep . Our starting point is to consider for each atom a set of K different observables {ô k (i)} K k=1 , some of which must be non-commuting in order to potentially detect entanglement. We then assume that data of the following form are available: where . . . = Tr(ρ . . . ) is an expectation value over many identically-prepared systems. In practice, as we illustrate below, the correlations C kl might be available for only a subset of pairs (k, l), something which can be accommodated within our approach. It is also central to our method to find the bound β defined as: where |ψ is an arbitrary internal state of atom i. Depending on the choice of observablesô k (i), the bound β might be computable analytically, but in general it may be found numerically, optimizing over the states of atom i.
In the examples studied in this work, the local observablesô k (i) are the same for all atoms (they are the components of the spin, or projectors onto spin states), and we simply parametrize |ψ = d α=1 c α |α with its complex coefficients c α in a fixed basis. K k=1 | ψ|ô k (i)|ψ | 2 is then a function of the coefficient c α (and is independent of i in the examples we study), that we optimized using standard numerical routines.
Our central result is the following inequality, valid for any separable state of the many-body system: where we introduced the notation [m ⊗ m] kl = m k m l , and P is any projector (namely a K × K symmetric matrix obeying P 2 = P ). The detailed proof of Eq. (4) is given in Appendix A. The optimal choice to most strongly violate Eq. (4) is to project onto the subspace corresponding to all negative eigenvalues of the matrix C − m ⊗ m. As a special case, Eq. (4) contains all generalized SSI of ref. [42] when choosing the spin in orthogonal directions x, y, z as local observables (see Appendix B for a detailed derivation). As we shall see, when applied to suitable states of a spinor gas, it allows us to find novel entanglement criteria beyond SSI.
Zeeman-sublevels population measurements. In order to illustrate the relevance of our new approach to probe entanglement in a spin-j atomic ensemble, we assume that the populationsN a,s of all Zeeman sublevels s = −j, . . . , j can be measured, for different quantization axes a (at least, for the purpose of entanglement detection, along two different orientations). This is typically achieved by first applying collective spin rotations to map a given orientation a onto the z axis, followed by a spatial splitting of the atoms in a magnetic-field gradient. The intensity of fluorescence imaging finally allows one to estimate the number of atoms in each Zeeman sublevel. For each orientation a, average populations N a,s and their correlations N a,sNa,s can be inferred in this manner. In the general framework presented above, the singleatom observables correspond to the projectors onto the spin statesn a,s (i) = |s a i s a | i (namely, the label k of previous section is now a pair of labels (a, s) denoting both the quantization axis and the spin state). Correspondingly, one reconstructs the average values m a,s = N a,s /N , and the diagonal blocks of the correlation matrix corresponding to a fixed orientation The projector P in Eq. (4) is then chosen in a block-diagonal form, projecting onto the negative eigenspace of each block, defined by the matrices C (a,s)(a,s ) − m a,s m a,s for the different quantization axes a.

Applications to a spin-1 BEC
We first illustrate our approach on a BEC of spin-1 atoms, as investigated recently in several experiments [9, 10, 18,33,45]. We consider N spin-1 atoms with contact isotropic interactions in the single-mode approximation, adiabatically prepared from the so-called polar state |0, N, 0 = ⊗ N i=1 |0 z i (with |N − , N 0 , N + the state written in the collective population basis along the z quantization axis). In an external magnetic field, the Hamiltonian reads (see Appendix C for details): whereĴ a = j s=−j sN a,s is the collective spin in direction a ∈ {x, y, z} andĴ 2 =Ĵ 2 x +Ĵ 2 y +Ĵ 2 z ;Q z = j s=−j s 2N z,s is the collective quadrupole operator; and {x, y, z} is an orthornormal basis of R 3 . Notice that the linear Zeeman term, proportional toĴ z , is omitted as it commutes withĴ 2 and acts trivially on the initial state: J z |0, N, 0 = 0. The initial state is also the ground state ofQ z ; and by varying the magnetic field intensity (q) one can adiabatically prepare the ground state ofĤ j=1 for either ferromagnetic interactions (c ≤ 0) or antiferromagnetic interactions (c ≥ 0). Notice that in the case of ferromagnetic interactions, two quantum phase transitions occur at q/|c| = ±4, so that adiabaticity is possible only for finite N . As the Zeeman sublevel populations are related by j s=−jN a, = N for all directions, m a,s and C (a,s)(a,s ) for all s, s are not all independent. As local observables in Eqs. (1) and (2) we used the spinŜ a (i) = j s=−j sn a,s (i) and the projector onto the Zeeman sublevel s = 0,n a,0 (i). We use these data as input to Eq. (4), with bound β j=1 = 3/2 in Eq. (3) (see Appendix F). Notice that for j = 1, we havê Q a =N a,+ +N a,− = N −N a,0 .
[10], where an almost perfect spin singlet is prepared, and all populations are measured; we also recovered Eq. (6) as the optimal violated EW. This EW is reminiscent of, yet generally more powerful than, the generalized SSI Var(Ĵ) ≥ N j (with j = 1 here) derived in ref. [42]. For a quantitative comparison of both criteria, in Fig. 1(a) we illustrate their detection power for N = 100 and states with Ĵ = 0 (vanishing mean spin, as is the case in the ground state ofĤ j=1 , Eq. (5)) in the plane (x = Ĵ 2 /N, y = N 2 0 /[N (N − 1)]), where violation occurs respectively for x − y < (N − 3)/[2(N − 1)] [Eq. (6)] and x < 1 (SSI). We plot the data points corresponding to the ground state ofĤ j=1 for c > 0 and both q > 0 and q < 0. More generally, not all combinations of Ĵ 2 and N 2 0 are physically allowed. As detailed in Appendix E (see also [38]), the physically-allowed region can be obtained by considering all quantum ground states of Hamiltonians of the form λ 1Ĵ 2 + λ 2N 2 0 for all values of λ 1 , λ 2 , and forming the convex envelope of the corresponding data points. The excluded (non-feasible) region is indicated in light grey in Fig. 1(a). Regarding the ground state ofĤ j=1 , three limiting cases are especially illustrative. For q → +∞, the ground state is the polar state |0, N, 0 = ⊗ N i=1 |0 z i , which is separable and therefore does not violate any EW. For q = 0, the ground state is a spin singlet, and violates both our witness Eq. (6) and the SSI Var(Ĵ) ≥ N j. In-between, we notice that our EW detects entanglement in the ground state ofĤ j=1 for a larger range of values of q than SSIs. Finally, for q → −∞, the ground state is a twin Fock state |N/2, 0, N/2 . This state does not violate any EW if only collective spin observables are measured (notice, though, that it does violate a generalized SSI for j = 1/2 particles, when considering only the two-dimensional subspace spanned by Zeeman sublevels s = ±1). Yet, the twin-Fock states robustly violates the EW of Eq. (6) (with Var(Ĵ) = N and N 2 0 = N (3N + 2)/4, so that the l.h.s is ∼ N/4, while the r.h.s is ∼ N/2). The polar, singlet and twin-Fock states are indicated on Fig. 1(a).

Ferromagnetic interactions.
For ferromagnetic interactions (c < 0) and q/|c| 1.2 (N = 100), we find violation of (see Appendix F): Maximal violation is found close to q = 0 where the state is a Dicke state, such thatĴ 2 |Dicke = N j(N j + 1)|Dicke (maximal total spin), and J z |Dicke = 0. Eq. (7) can be compared to another SSI of ref. [42] violated by the Dicke state, Fig. 1(b), we compare their power in detecting entanglement in the ground state ofĤ j=1 . In particular, for q |c|, the ground state ofĤ j=1 is very close to the twin-Fock state |N/2, 0, N/2 , which robustly violates Eq. (7): the l.h.s is equal to N (N + 1)/2 ∼ N 2 /2, while on the r.h.s N 2 0 = N (3N + 2)/4 ∼ 4N 2 /4. To conclude on the spin-1 BEC, we find that collective population measurements in orthogonal Zeeman sublevels can reveal entanglement beyond collective spin observables, and our method can be used to infer novel EWs in an optimal and data-agnostic way from these measurements. These new EWs are generically more powerful than SSIs (which can be recovered as a special case in our approach, see Appendix B). In particular, we exhibit two EWs [Eqs. (6) and (7)] robustly violated by the twin-Fock state |N/2, 0, N/2 , stabilized as the ground state of H j=1 in the limit q |c|; yet no generalized SSI for spin-1 ensembles is violated by the twin-Fock state.

Applications to a spin-2 BEC
The Hamiltonian of a BEC of j = 2 atoms interacting via pairwise isotropic interactions reads (see Appendix C): In contrast to the spin-1 case [Eq. (5)], the Hamiltonian contains the additional pairing termΘ † 2Θ 2 whereΘ 2 = j s=−j (−1) sâ sâ−s annihilates timereversal pairs ±s, withâ s annihilating an atom in the spin state |s z . Notice that [Ĵ 2 ,Θ † 2Θ 2 ] = 0. We explored the ground state ofĤ j=2 to find novel violated EW. Using as local operatorsŜ a (i) (the spin) andn a,0 (i) (the projector onto the s = 0 spin state) for a ∈ {x, y, z}, the bound in Eq.
(3) is numerically found as β 2 ≈ 4.3. For this choice of local operators, we found that generalizations of Eqs. (6) and (7) are violated. These EW can be further generalized to arbitrary integer spins, and read (see Appendix F): and: where 1 = (1, 1, 1) and β j is the corresponding bound in Eq.
(3). As for the spin-1 case, one can compare the entanglement detection power of these witnesses with the related SSIs; this comparison is presented in Fig. 2 for Eq. (9), while for Eq. (10) we notice that the improvement with respect to the SSI is mild, and not systematic. Interestingly, as in the case of spin-1, the twin-Fock state |N/2, 0, 0, 0, N/2 (obtained in the limit q → −∞) violates Eq. (9), while it does not violate any generalized SSI for spin-2 collective spin observables. As a final exploration of the entanglement patterns in the ground state of a spin-2 BEC, we have focused on many-body singlets ( Ĵ 2 = 0), stabilized at q = 0 for antiferromagnetic interactions (c > 0). In the spin-1 case, there is a unique such state for a BEC, obtained as a "pair condensate" (Θ † 2 ) N/2 |vac. . This pair condensate is also stabilized in the spin-2 case for p < 0. Yet, for p > 0, the many-body singlet is instead a "triplet condensate" is the singlet-trio creation operator. Both states are entangled and violate Eq. (9). However, considering as local measurements {n a,+2 (i),n a,+1 (i),n a,−1 (i),n a,−2 (i)} a∈{x,y,z} , we found the following EW violated only by the pair condensate and not by the triplet condensate (see Appendix G):

Conclusions
We have presented a new method to infer violated entanglement witnesses from the collective measurement of Zeeman-sublevel populations in spinor gases of arbitrary spin-j atoms.
Our method recovers all known generalized spin squeezing inequalities [42] as a special case, in a data-agnostic way. But when considering appropriate correlations between Zeeman-sublevel populations, one can also construct novel entanglement criteria which have no analog if only the collective spin is measured. We have illustrated our method on the ground states of spin-1 and spin-2 BECs, interacting with contact twobody isotropic interactions in an external magnetic field. In particular, for both j = 1 and j = 2 we showed that the twin Fock state robustly violates an entanglement witness involving all populations in three orthogonal directions, while entanglement cannot be detected if only the collective spin is measured. For j = 2, we presented an entanglement witness distinguishing between two forms of many-body singlets, namely a pair condensate and a triplet condensate. Overall, our new approach appears especially suited to probe entanglement among many particles in high-spin spinor gases, for which entanglement detection has remained a considerable challenge. We would like to emphasized that, although we have focused on collective measurements which are invariant under the permutations of the atoms, our approach is not limited to this situation. It could also be used to probe entanglement in spatiallystructured systems by using different local observablesô k (i) for each atom; a straightforward example is to introduce local phases e iφ k (i) , leaving the bound β invariant in Eq.
(3). As pointed out recently [12], this offers considerable flexibility to probe entanglement through components of the structure factors. Finally, a limitation of our current implementation is that there is no optimal and systematic way of incorporating all available first-and second-order moments of population measurements into our main inequality (4). Indeed, while adding more measurement operators a priori improves the entanglement detection power, this leads to an increase of the β bound in Eq.

A Derivation of the main entanglement witness inequality
In this section, we derive our central result, Eq. (4). In order to detect entanglement from the available data [Eqs.
(1) and (2)], our starting point to is to assume that the latter are explained by a separable state, namely a quantum state of the form: and to arrive at a contradiction. We introduce the notations: The key insight of our approach is that, if the state is separable, we may interpret the quantities m k (λ, i) as local classical variables, whose collective fluctuations must reproduce the observed fluctuations of collective observables. The constraints that such local classical variables must obey give rise to our central inequality [Eq. (4)]. Most importantly, the constraints of Eq.
(3) can be expressed as: for all λ. With the above notations, the one-body quantum data [Eq.
(1)] read: We then introduce the vector notation m := (m k ), and the outer product notation With these notations, the two-body quantum data [Eq.
(2)] read: We now introduce the fluctuations of the m k (λ, i) variables: The notation M 0 means the the matrix M is positive semidefinite (PSD), that is, for any vector u, u † M u ≥ 0. This PSD property is straightforward to verify, since A and a are both of the form The matrix C is then: From Eq. (20) and Eq. (15)], we have: We then introduce a projector P , and write:

= Tr[P (m * ⊗ m + A)] + Tr[(1 − P )(C + a/(N − 1))] (24)
where we used Eq. (21). From the positivity of A and P , we have Tr(P A) ≥ 0. Similarly, from the positivity of a and 1 − P (which is also a projector), we have Tr[(1 − P )a] ≥ 0. Hence: Introducing this inequality into Eq. (22), and reorganizing the terms, this leads to our central result: which is valid for any projector P whenever the state is separable [Eq. (12)]. This central inequality is a generalization of the results contained of Vitagliano et al. [42]. In order to detect entanglement at best, one has to find the optimal projector P such that inequality (26) is violated -which would invalidate our starting assumption, namely that the state is separable. This optimization is straightforward: in order to minimize L over the projector P , one has to choose P as the projector onto the subspace of negative eigenvalues of C − m * ⊗ m.
In ref. [42], a family of SSI was derived that can be obtained as a special case of our approach. These inequalities are valid for all fully-separable states of N spin-j atoms, and involve first and (modified) second moments of collective spins in three orthogonal directions:Ĵ x ,Ĵ y ,Ĵ z . We write here these inequalities (Eq. (9) of ref. [42]) for completeness: The labels k, l, m in the last two inequalities may take all permutations of x, y, z. As we shall show, these inequalities are obtained from our main inequality Eq. (4) (or equivalently Eq. (26) in Appendix A), using as local observablesô k (i) the three spin observablesŜ x (i),Ŝ y (i),Ŝ z (i). We then have: Furthermore, the β bound is β = j 2 , since for any state of a spin-j atom, one has S x (i) 2 + S y (i) 2 + S z (i) 2 ≤ j 2 , with equality if the atomic spin is fully polarized along some direction. Hence we have: where we used that [S x (i)] 2 + [S y (i)] 2 + [S z (i)] 2 = j(j + 1).
with E any subset of {x, y, z}. This last inequality can be written more explicitly as: Choosing E = ∅ directly gives Eq. (27).
To obtain Eqs. (29) and (30), we rewrite Eq. (35) as: Taking then e.g. E = {x}, we find: namely Eq. This concludes the proof that the generalized SSI of ref. [42] are incorporated in the framework presented in this paper when choosing the spin components as local observables.

C The spinor Bose gas model
In this section, we present some background information, based on Ref. [21], on the BEC Hamiltonians used in the main text [Eqs. (5) and (8)]. We consider an ensemble of spin-j bosonic atoms with contact, two-body, spin-conserving interactions in R 3 . In the single-mode approximation, its Hamiltonian reads: where g f is the scattering amplitude of the spin-2f channel. The matrix element of the projection of a composite spin-2f onto two spin-j atoms is C j⊗j=2f s 1 s 2 ,s 3 s 4 = j, s 1 | j ⊗ j, s 2 | j P j⊗j=2f |j, s 3 j ⊗ |j, s 4 j , with P j⊗j=2f = 2f S=−2f |2f, S j⊗j 2f, S| j⊗j , where S is the spin projection of the composite spin-2f . The operatorsâ s (â † s ) annihilate (create) a bosonic mode with spin projection s ∈ {−j, −j + 1, .., j}. The total number of atomsN = j s=−jN s , whereN s =â † sâ s , is fixed to N . In the following we will express the Hamiltonian (39) in terms of collective operators. This is done via the identification:Ô whereÔ = N i=1ô (i) is the collective operator corresponding to the local observableô acting on a single-atom spin-j degree of freedom. Note that in particular the population operatorN s is a collective operator with s L |n s |s R := s L |s s|s R = δ s L ,s δ s,s R . Spin 1. For spin-1, the spin projection onto the cartesian basis {x, y, z} reads: where i is the imaginary unit, i 2 = −1. With this expressions, we can compute the total spin, using the bosonic commutation relations and N z,+1 + N z,0 + N z,−1 = N we have: Alternatively, by identifying the pair annihilation operatorΘ 2 =â 2 0 − 2â +1â−1 , we can writeĴ 2 = N 2 −Θ †Θ . As first observation, we notice that unlike the spin-1/2 case,Ĵ 2 is not diagonal in the occupation basis. Direct computation shows that the Hamiltonian (39) is equivalent to: where c 1 = (g 1 − g 0 )/3.
Quadratic Zeeman shift. In this work we consider the gas is immersed in a uniform magnetic field of intensity B along the z direction. In this case, the total Hamiltonian becomesĤ We introduce a scaling to make the Hamiltonian extensive in N . We consider an experiment where all atoms are initially prepared in the state ⊗ N i=1 |0 z , so that we may remove the linear Zeeman term, proportional toĴ z (asĴ z is a symmetry ofĤ, and acts trivially on the state).

D Computation of the quantum data
In this appendix, we provide technical details on the computation of the quantum data used in the illustration of our method. Specifically, one has to compute one and two body Zeeman sublevel correlations of the form: N a,s , N a,sNa,s , etc. , in the ground state of relevant Hamiltonians of spin-j ensembles.
Diagonalization of the Hamiltonians. The Hamiltonians considered in this paper are polynomials of collective operatorsÔ k = N i=1ôk (i) with the spin projectionĴ z a symmetry. Since operators acting in different atoms commute, we notice that the collective operators {Ô k } obey the same algebra as their their local counterparts, {ô k (i)}. Accordingly, the problem reduces to compute the matrices of the irreducible representations (irreps) of the algebra, describing the collective operators. The local operators can be described by d × d Hermitian matrices, where d = 2j + 1. Let {λ α } d 2 −1 α=1 be a basis of d × d Hermitian traceless matrices fulfilling the orthogonality condition Tr[λ αλβ ] = δ αβ , for instance the so-called generalized Gell-Mann matrices. The commutation relations: Consequently, for any other operatorô, the map I distributes linearly where D is the dimension of the representation, and the identity element I is added in order to account for the trace ofô. The multiplicity N is introduced to consistently transform the identity, (49) is an alternative of the second-quantization form (40) providing the advantage that higher moments are computed just by multiplying matrices, as specified in the following.

E Computation of the physically feasible region of a set of average values
In this section, we explain how to compute the feasible region showed in the figures of the main text. Formally, given two observablesX andŶ , acting on the Hilbert space for N qudits, our goal is find the possible combinations of "feasible" mean values (x = X , y = Ŷ ), that is, combinations of mean values which can be obtained in some (generally mixed) quantum state. Equivalently, one may look for the extremal values of Ŷ , while keeping fixed the mean value x = X . By making statistical mixtures of the corresponding quantum states, the full convex hull of these extremal points is also feasible. This problem can be solved introducing a Lagrange multiplier λ x , and finding the ground state of Hamiltonians of the form:L Adjusting the multiplier λ x ensures that the constraint X = x is met in the ground state. More generally, the extremal points (x, y) are found as the expectation values ofX,Ŷ in the ground state of:L varying the parameters λ x , λ y between −∞ and +∞. By forming statistical mixtures of these ground states, the full convex hull of the extremal points are also feasible, which leads to the construction showed in the figures of the main text. In order to find the ground state ofL, the same technique outlined above for the Bose gas Hamiltonian can be applied since in the cases of interest here,L is a polynomial of collective operators.
F Derivation of the entanglement witness inequalities (6), (7) and their generalization to arbitrary integer spins Eqs. (9), (10) These inequalities are found by choosing, as local measurements, the spin observablesŜ a (i) = j s=−j sn a,s (i), as well as the projectorsn a,0 (s) onto the s = 0 magnetic sublevels, along the three orthogonal directions a ∈ {x, y, z}. For completeness, we rewrite here the EW inequalities (9) and (10) of the main text: and: where 1 = (1, 1, 1) and β j is the corresponding bound in Eq.
Derivation of Eq. (52). For Eq. (52) (tailored in particular to a spin-j many-body singlet), the optimal projection P is P x = P y = P z = 1 0 0 0 := P j in the local basis described before, leading to Inserting this expression, together with Eq. (56), into our central inequality (4), we obtain the announced result Eq. (52).
Derivation of Eq. (53). For Eq. (53) (tailored in particular to the Dicke states J 2 = N j(N j + 1), J z = 0), the optimal projector is P x = P y = 0, P z = P j , from which we find: Inserting this expression, together with Eq. (56), into our central inequality (4), we obtain the announced result Eq. (53).
G Derivation of spin-2 EW, inequality (11) and generalization to arbitrary integer spins In this appendix we show how to derive the entanglement witness presented in Eq. (11), and violated by the pair condensate, but not by the triplet condensate. In the present case, we use as local observables the projectors onto all magnetic sublevels, except the s = 0 sublevel: {n a,+2 (i),n a,+1 (i),n a,−1 (i),n a,−2 (i)} a∈{x,y,z} . This choice leads to the Tr(C) term: We introduce the notationσ a,s =N a,s −N a,−s , the population imbalance between opposed levels along direction a. Together with the notation 1 = (1, 1, 1), we can rewrite these equations as: namely, Eq. (11).

Generalization to arbitrary even-level systems.
A possible generalization to arbitrary integer spin is obtained by taking the projector to j s=1 (n a,s −n a,−s ) for a ∈ {x, y, z}. The corresponding bound β for local measurements {{n a,s } j s=−j =0 } a∈{x,y,z} obtained numerically is summarized in the following table for several spin-j.  Notice that for j = 1, we recover Eq. (6).

H Practical implementation of the algorithm
In this appendix, we review step by step the implementation of the presented algorithm in order to reveal quantum entanglement from experimental data. The procedure consists of two parts. (i) First, we infer the required mean values from some hypothetical experiment. Then (ii), we apply our central result Ineq. (4) to (hopefully) detect entanglement. Note that our approach is system-agnostic, in the sense that it does not depend on features of the system not captured by the correlations considered. However, for the sake of concreteness, we will illustrate the implementation of the method with a concrete example. Consistently with the applications presented in the main text, we will consider a BEC of spin-1 atoms.

H.1 Setting up the scenario
The first step is to specify the partitioning of the considered quantum system to define entangled states. Specifically, here we consider an ensemble of N three-level atoms; and we define a state as entangled if it cannot be decomposed as a mixture of product states over single-atom wave-functions. The next step is to specify the single-atom observables {{ô k (i)} K k=1 } N i=1 . Entanglement will be detected from correlations among them, so it must contain at least a pair of non-commuting observables. Furthermore, we consider the same set of local measurements for each atom. This simplification will allow us to infer the correlations in a scalable way from moments of collective operators, which is also the natural way in which measurements are performed in many cold-atom experiments. Concretely, in the present paper we consider two classes of local operators: (i) spin observables; and (ii) projectors to Zeeman sublevels. Here we work with {{n a,+1 (i),n a,−1 (i)} K a={x,y,z} } N i=1 , wheren a,s (i) is the projector to the hyperfine level s along direction a for atom i as defined in the main text and {x, y, z} an orthonormal basis of R 3 . Using a certain data to be specified in the next subsection, this setting leads to the witness Ineq. (6) of the main text, which we originally derived from a different set of local observables, namely for a combination of spin and populations to the zero state {{ŝ a (i),n a,0 (i)} K a={x,y,z} } N i=1 .