Guillem Müller-Rigat

Maciej Lewenstein

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 fluctuations. 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 finer 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 specific 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 dlevel 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 [START_REF] Gilmore | Quantum-enhanced sensing of displacements and electric fields with two-dimensional trapped-ion crystals[END_REF][START_REF] Pezzè | Quantum metrology with nonclassical states of atomic ensembles[END_REF], or the simulation of quantum many-body problems [START_REF] Georgescu | Quantum simulation[END_REF].

Guillem Müller-Rigat: guillem.muller@icfo.eu Irénée Frérot: irenee.frerot@neel.cnrs.fr 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 [START_REF] Kitagawa | Squeezed spin states[END_REF][START_REF] Ma | Quantum spin squeezing[END_REF][START_REF] Pezzè | Quantum metrology with nonclassical states of atomic ensembles[END_REF][START_REF] Sørensen | Many-particle entanglement with bose-einstein condensates[END_REF][START_REF] Sørensen | Entanglement and extreme spin squeezing[END_REF][START_REF] Tóth | Spin squeezing and entanglement[END_REF][START_REF] Vitagliano | Spin squeezing inequalities for arbitrary spin[END_REF][START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF][START_REF] Wineland | Squeezed atomic states and projection noise in spectroscopy[END_REF]. SSI and their generalizations [START_REF] Gessner | Metrological nonlinear squeezing parameter[END_REF][START_REF] Tóth | Spin squeezing and entanglement[END_REF][START_REF] Vitagliano | Spin squeezing inequalities for arbitrary spin[END_REF][START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF] have found an impressively broad range of successful applications, from detecting entanglement in cold and hot atomic ensembles [START_REF] Kong | Measurement-induced, spatially-extended entanglement in a hot, strongly-interacting atomic system[END_REF][START_REF] Pezzè | Quantum metrology with nonclassical states of atomic ensembles[END_REF] to unveiling many-body Bell non-locality [START_REF] Engelsen | Bell Correlations in Spin-Squeezed States of 500 000 Atoms[END_REF][START_REF] Müller-Rigat | Inferring nonlinear many-body bell inequalities from average two-body correlations: Systematic approach for arbitrary spinj ensembles[END_REF][START_REF] Piga | Bell correlations at ising quantum critical points[END_REF][START_REF] Schmied | Bell correlations in a Bose-Einstein condensate[END_REF], from probing quantum phase transitions [START_REF] Frérot | Quantum critical metrology[END_REF] and quantum quenches [START_REF] Comparin | Multipartite entangled states in dipolar quantum simulators[END_REF][START_REF] Comparin | Robust spin squeezing from the tower of states of u(1)-symmetric spin hamiltonians[END_REF] to applications in quantum-enhanced metrology [START_REF] Pezzè | Quantum metrology with nonclassical states of atomic ensembles[END_REF]. 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 [START_REF] Tóth | Spin squeezing and entanglement[END_REF]. 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 [START_REF] Duan | Quantum entanglement in spinor boseeinstein condensates[END_REF][START_REF] Hamley | Spin-nematic squeezed vacuum in a quantum gas[END_REF][START_REF] Kunkel | Spatially distributed multipartite entanglement enables epr steering of atomic clouds[END_REF][START_REF] Kunkel | Simultaneous readout of noncommuting collective spin observables beyond the standard quantum limit[END_REF][START_REF] Qu | Probing spin correlations in a bose-einstein condensate near the single-atom level[END_REF]. 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 [START_REF] Rohit | Stochastic dynamics of a few sodium atoms in presence of a cold potassium cloud[END_REF][START_REF] Kasper | Universal quantum computation and quantum error correction with ultracold atomic mixtures[END_REF], or effective qudits ensembles [START_REF] Ringbauer | A universal qudit quantum processor with trapped ions[END_REF]. 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. [START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF] 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. [START_REF] Vitagliano | Spin squeezing inequalities for arbitrary spin[END_REF], 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 [START_REF] Tóth | Spin squeezing and entanglement[END_REF][START_REF] Vitagliano | Spin squeezing inequalities for arbitrary spin[END_REF][START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF], 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 [START_REF] Evrard | From many-body oscillations to thermalization in an isolated spinor gas[END_REF][START_REF] Evrard | Observation of fragmentation of a spinor bose-einstein condensate[END_REF][START_REF] Hamley | Spin-nematic squeezed vacuum in a quantum gas[END_REF][START_REF] Qu | Probing spin correlations in a bose-einstein condensate near the single-atom level[END_REF][START_REF] Zou | Beating the classical precision limit with spin-1 dicke states of more than 10,000 atoms[END_REF]. 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:

m k = 1 N N i=1 ôk (i) , (1)
C kl = 1 N (N -1) i =j ôk (i)ô l (j) , (2)
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:

β = 1 N N i=1 max |ψ K k=1 | ψ|ô k (i)|ψ | 2 , ( 3
)
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:

N Tr[P (C -m ⊗ m)] ≥ Tr(C) -β , (4)
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. [START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF] 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 populations Na,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 Na,s and their correlations Na,s Na,s can be inferred in this manner. In the general framework presented above, the singleatom observables correspond to the projectors onto the spin states na,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 = Na,s /N , and the diagonal blocks of the correlation matrix corresponding to a fixed orientation

C (a,s)(a,s ) = Na,s [ Na,s -δ s,s ] /[N (N -1)].
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 [START_REF] Evrard | From many-body oscillations to thermalization in an isolated spinor gas[END_REF][START_REF] Evrard | Observation of fragmentation of a spinor bose-einstein condensate[END_REF][START_REF] Hamley | Spin-nematic squeezed vacuum in a quantum gas[END_REF][START_REF] Qu | Probing spin correlations in a bose-einstein condensate near the single-atom level[END_REF][START_REF] Zou | Beating the classical precision limit with spin-1 dicke states of more than 10,000 atoms[END_REF]. 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):

Ĥj=1 = c N Ĵ2 + q Qz , ( 5
)
where Ĵa = j s=-j s Na,s is the collective spin in direction a ∈ {x, y, z} and Ĵ2 = Ĵ2

x + Ĵ2 y + Ĵ2 z ; Qz = j s=-j s 2 Nz,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: Ĵz |0, N, 0 = 0. The initial state is also the ground state of Qz ; 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=-j Na, = 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, na,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 have Qa = Na,+ + Na,-= N -Na,0 . Antiferromagnetic interactions. For antiferromagnetic interactions (c > 0) and q/c 4.4 (N = 100), we find violation of Eq. (4), with an entanglement witness (EW) which reads (see Appendix F):

Var( Ĵ) - 1 N -1 N2 0 ≥ N (N -3) 2(N -1) , ( 6
)
with the notations Var( Ĵ) = Var( Ĵx ) + Var( Ĵy ) + Var( Ĵz ), Var( Ô) = Ô2 -Ô 2 and N0 = ( Nx,0 , Ny,0 , Nz,0 ). The EW of Eq. ( 6) is maximally violated at q = 0, where the ground state of Eq. ( 5) is a spin singlet (Var( Ĵ) = 0). In particular, we also applied our method using as input the experimental data of ref. [START_REF] Evrard | Observation of fragmentation of a spinor bose-einstein condensate[END_REF], 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. [START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF]. 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 6)] (solid line) with the SSI Var( Ĵ) ≥ N j (dashed line) for spin-1 states (N = 100 atoms). The light grey area are data points which cannot be achieved (see text and Appendix E). The dotted line indicates the ground state of Eq. ( 5) for both antiferromagnetic interactions (c = 1), varying the quadratic Zeeman term (-10 ≤ q ≤ 10). Three representative states are also indicated: the polar state |0, N, 0 (q → ∞) which is separable, the spin singlet s.t. Ĵ2 |Ψ = 0 (q = 0), and the twin-Fock state |N/2, 0, N/2 (q → -∞). (b) Comparison of the population-based witness [Eq. [START_REF] Duan | Quantum entanglement in spinor boseeinstein condensates[END_REF]] with a related SSI (see text) for spin-1 states. Graphical conventions as in panel (a). Here, the ground state is computed for ferromagnetic interactions (c = -1). At q = 0, we obtain the Dicke state

Ĵ2 |Dicke = N (N + 1)|Dicke , Ĵz |Ψ = 0. plane (x = Ĵ2 /N, y = N2 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 N2 0 are physically allowed. As detailed in Appendix E (see also [START_REF] Sørensen | Entanglement and extreme spin squeezing[END_REF]), the physically-allowed region can be obtained by considering all quantum ground states of Hamiltonians of the form

λ 1 Ĵ2 + λ 2 N2
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 N2 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):

(N -1)Var( Ĵz ) -N Qz -Ĵ2 x + Ĵ2 y + 3N (N + 1) 2 ≥ N2 0 . ( 7
)
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 Ĵz |Dicke = 0. Eq. ( 7) can be compared to another SSI of ref. [START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF] violated by the Dicke state, namely (N -1)Var( Ĵz ) -N Qz -Ĵ2

x + Ĵ2 y + N j(N j + 1) ≥ 0. In 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 N2 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. [START_REF] Comparin | Robust spin squeezing from the tower of states of u(1)-symmetric spin hamiltonians[END_REF] and [START_REF] Duan | Quantum entanglement in spinor boseeinstein condensates[END_REF]] robustly violated by the twin-Fock state |N/2, 0, N/2 , stabilized as the ground state of Ĥ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):

Ĥj=2 = c N Ĵ2 + p N Θ † 2 Θ2 + q Qz . ( 8
)
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) and na,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:

Var( Ĵ) - 1 N -1 N0 ( N0 -1) ≥ N [j(j + 1) -β j ] , (9)
(N -1)Var( Ĵz ) -N Qz -Ĵ2 x + Ĵ2 y +N [(N -1)β j + j(j + 1)] ≥ N0 ( N0 -1) , ( 10
)
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. [START_REF] Evrard | Observation of fragmentation of a spinor bose-einstein condensate[END_REF] 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" ( Θ †

3 ) N/3 |vac. , where

Θ † 3 = -9[(a † -1 ) 2 a † 2 +(a † 1 ) 2 a † -2 ]+ √ 6a † 0 [-(a † 0 ) 2 +3a † 1 a † -1 + 6a † 2 a † -2 ]
is the singlet-trio creation operator. Both states are entangled and violate Eq. (9). However, considering as local measurements {n a,+2 (i), na,+1 (i), na,-1 (i), na,-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):

σ2 2 + σ2 1 - N -1 N ( σ2 2 + σ1 2 ) - 2 N N2 +2 + N2 +1 + N2 -1 + N2 -2 +(1 - 2 N ) N0 • 1 + 3 ≥ 0 , ( 11
)
where we denoted by σa,s = Na,s -Na,-s the population imbalance between opposed levels along direction a. Notice that {σ x,s , σy,s , σz,s } generate an su(2) subalgebra of su( 5). This result illustrates the difference in entanglement properties of two different spin singlets, which by definition cannot be distinguished if only collective spin observables are measured.

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 [START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF] 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 [START_REF] Frérot | Unveiling quantum entanglement in many-body systems from partial information[END_REF], 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. (3), and can make Eq. ( 4) harder to violate. The optimization of the choice of measurement operators is left open to future studies.

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:

ρsep = λ p λ ⊗ N i=1 ρλ (i) , ( 12
)
and to arrive at a contradiction. We introduce the notations:

m k (λ, i) := Tr[ρ λ (i)ô k (i)] ( 13
)
m k (λ) := 1 N N i=1 m k (λ, i) . ( 14
)
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:

1 N N i=1 k [m k (λ, i)] 2 ≤ β (15)
for all λ. With the above notations, the one-body quantum data [Eq. (1)] read:

m k = 1 N N i=1 ôk (i) = λ p λ m k (λ) (16)
We then introduce the vector notation m := (m k ), and the outer product notation

[u ⊗ v] ab = u a v b .
With these notations, the two-body quantum data [Eq. (2)] read:

C = 1 N (N -1) i =j λ p λ m(λ, i) * ⊗ m(λ, j) . ( 17
)
We now introduce the fluctuations of the m k (λ, i) variables:

A := λ p λ [m(λ) -m] * ⊗ [m(λ) -m] 0 (18)
a := λ p λ 1 N i [m(λ, i) -m(λ)] * ⊗ [m(λ, i) -m(λ)] 0 (19)
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 x w x vv † with w x ≥ 0 and v a one-dimensional vector. The matrix C is then:

C = N N -1 [m * ⊗ m + A] - 1 N (N -1) λ p λ i m(λ, i) * ⊗ m(λ, i) (20) = m * ⊗ m + A - 1 N -1 a . ( 21
)
From Eq. (20) and Eq. ( 15)], we have:

Tr(C) ≥ N N -1 Tr[m * ⊗ m + A] - β N -1 . ( 22
)
We then introduce a projector P , and write:

Tr[m * ⊗ m + A] = Tr[P (m * ⊗ m + A)] + Tr[(1 -P )(m * ⊗ m + A)] (23) = Tr[P (m * ⊗ m + A)] + Tr[(1 -P )(C + a/(N -1))] ( 24
)
where we used Eq. [START_REF] Kawaguchi | Spinor bose-einstein condensates[END_REF]. 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:

Tr[m * ⊗ m + A] ≥ Tr[P (m * ⊗ m)] + Tr[(1 -P )C] ( 25
)
Introducing this inequality into Eq. ( 22), and reorganizing the terms, this leads to our central result:

L(P ) := N Tr[P (C -m * ⊗ m)] ≥ Tr(C) -β , ( 26
)
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. [START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF]. 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.

B Recovering the generalized spin-squeezing inequalities of ref. [START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF].

In ref. [START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF], 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. [START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF]) for completeness:

Ĵ2 x + Ĵ2 y + Ĵ2 z ≤ N j(N j + 1) (27)
Var( Ĵx ) + Var( Ĵy ) + Var( Ĵz ) ≥ N j (28)

Ĵ2 l - N i=1 [ Ŝl (i)] 2 + Ĵ2 m - N i=1 [ Ŝm (i)] 2 ≤ (N -1) Var( Ĵk ) - N i=1 [ Ŝk (i)] 2 +N (N -1)j 2 (29) Ĵ2 m - N i=1 [ Ŝm (i)] 2 -N (N -1)j 2 ≤ (N -1) Var( Ĵk ) - N i=1 [ Ŝk (i)] 2 +Var( Ĵl ) - N i=1 [ Ŝl (i)] 2 (30)
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:

m k = Ĵk /N ( 31
)
C kl = 1 N (N -1) Ĵk Ĵl - i Ŝk (i) Ŝl (i) . ( 32
)
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:

Tr(C) -β = 1 N (N -1) ( J 2 x + J 2 y + J 2 z ) - 1 N -1 j(j + 1) -j 2 (33)
where we used that [S x (i)] 2 + [S y (i)] 2 + [S z (i)] 2 = j(j + 1). As we will show explicitly, the 8 SSI [Eqs. [START_REF] Lucioni | Dysprosium dipolar bose-einstein condensate with broad feshbach resonances[END_REF] to [START_REF] Müller-Rigat | Inferring nonlinear many-body bell inequalities from average two-body correlations: Systematic approach for arbitrary spinj ensembles[END_REF]] are obtained by choosing the projector P in Eq. ( 26) diagonal, leading to inequalities of the form:

N k∈E [C kk -m 2 k ] ≥ Tr(C) -j 2 (34)
with E any subset of {x, y, z}. This last inequality can be written more explicitly as:

k∈E N Ĵ2 k -(N -1) Ĵk 2 -N i [S k (i)] 2 ≥ J 2 x + J 2 y + J 2 z -N j(N j + 1) . ( 35
)
Choosing E = ∅ directly gives Eq. [START_REF] Lucioni | Dysprosium dipolar bose-einstein condensate with broad feshbach resonances[END_REF].

Choosing E = {x, y, z} gives:

(N -1) Var( Ĵx ) + Var( Ĵy ) + Var( Ĵz ) ≥ N j(j + 1) -N j(N j + 1) , [START_REF] Sørensen | Many-particle entanglement with bose-einstein condensates[END_REF] namely Eq. [START_REF] Luo | Deterministic entanglement generation from driving through quantum phase transitions[END_REF].

To obtain Eqs. ( 29) and (30), we rewrite Eq. ( 35) as:

k∈E N Ĵ2 k -(N -1) Ĵk 2 -N i [S k (i)] 2 ≥ J 2 x + J 2 y + J 2 z - i [ Ŝx (i)] 2 + [ Ŝy (i)] 2 + [ Ŝz (i)] 2 -N (N -1)j 2 . ( 37
)
Taking then e.g. E = {x}, we find:

(N -1) var( Ĵx ) - i [ Ŝx (i)] 2 ≥ J 2 y + J 2 z - i [ Ŝy (i)] 2 + [ Ŝz (i)] 2 -N (N -1)j 2 , ( 38
)
namely Eq. ( 29) for (k, l, m) = (x, y, z), and similarly for E = {y} or E = {z}. Finally, Eq. ( 30) is obtained by choosing E = {x, y}, E = {x, z} or E = {y, z}. This concludes the proof that the generalized SSI of ref. [START_REF] Vitagliano | Spin squeezing and entanglement for an arbitrary spin[END_REF] 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. [START_REF] Kawaguchi | Spinor bose-einstein condensates[END_REF], on the BEC Hamiltonians used in the main text [Eqs. [START_REF] Comparin | Multipartite entangled states in dipolar quantum simulators[END_REF] and [START_REF] Engelsen | Bell Correlations in Spin-Squeezed States of 500 000 Atoms[END_REF]]. 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:

Ĥ0 ({g f } j f =0 ) = j s 1 s 2 ,s 3 s 4 =-j g f C j⊗j=2f s 1 s 2 ,s 3 s 4 â † s 1 â † s 2 âs 3 âs 4 , (39)
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 atoms N = j s=-j Ns , where Ns = â † 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:

Ô = j s L ,s R =-j s L | ô |s R â † s L âs R , ( 40
)
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 operator Ns is a collective operator with

s L | ns |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:

Ĵx = 1 √ 2 â+1 â0 â-1    0 1 0 1 0 1 0 1 0       â+1 â0 â-1    (41)
Ĵy = 1 √ 2 â+1 â0 â-1    0 -i 0 i 0 -i 0 i 0       â+1 â0 â-1    (42)
Ĵz = â+1 â0 â-1

   1 0 0 0 0 0 0 0 -1       â+1 â0 â-1    , ( 43
)
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:

Ĵ2 : = Ĵ2 x + Ĵ2 y + Ĵ2 z = = Ĵ2 z + N + N0 + 2 N0 (N -N0 ) + 2((â † 0 ) 2 â+1 â-1 + h.c.) (44)
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:

Ĥ0 (c 1 ) = c 1 Ĵ2 , (45)
where c 1 = (g 1 -g 0 )/3.

Spin 2.

For spin-2, there is a second relevant interaction channel. In this case, the Hamiltonian reads:

Ĥ0 (c 2 , p) = c 2 Ĵ2 +p[ N0 ( N0 -1)+4( N1 N-1 + N-2 N2 )+2( T + 0,-2 T + 0,2 -T + 0,-1 T + 0,1 -2 T + 1,2 T + -1,-2 +(+ ↔ -))] ,
(46) where c 2 = (g 4 -g 2 )/7, p = (7g 0 -10g 2 + 3g 4 )/35 and we introduced the collective operators T + s,s = â † s âs + h.c., Ts,s = i(â † s âsh.c.). Using the pair anhiliation operator for spin-2, Θ2 = â2 0 -2â 1 â-1 + 2â +2 â-2 , the term coupling p can be interpreted as a pair amplitude Θ † 2 Θ2 . This term commutes with Ĵ2 , although we do not have written it explicitly in a rotationally invariant way using cartesian coordinates.

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 Ĥ = Ĥ0 /N + q Qz , where q ∝ B 2 and Qz = N i=1 Ĵ2 z . 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:

Na,s , Na,s Na,s , etc. , (47)

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. α=1 . The operator Λα is the collective counterpart of λα , that is Λα = N i=1 λα (i). Consequently, for any other operator ô, the map I distributes linearly

I[ô] = N Tr[ô] I D + d 2 -1 α=1 Tr[ λα ô] Λα := Ô , ( 49
)
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,

I[1] = N i=1 1 = N .
The expression (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.

Gathering the correlations. The Hamiltonian is diagonalized using the previous matrix representation, allowing one to find the ground state |GS . Then, the computation of moments reduces to matrix multiplication matrices and computing expectation values in the ground state: Ôa = GS| I[ô a ] |GS , Ôa Ôb = GS| I[ô a ]I[ô b ] |GS , etc, for any collective operators Ôa . In this work in particular we probe the GS with magnetic-sublevel populations along a direction a, { Na,s } j s=-j . These are generated from the single-atom observable na,s = j s =s =-j (s -Ŝa )/(s -s ), where Ŝa is the atomic spin projection along direction a ∈ {x, y, z}, i.e Ŝa := a • Ŝ.

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 observables X 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(λ x ) = ± Ŷ -λ x ( X -x1) . ( 50
)
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 of X, Ŷ in the ground state 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

Tr[N P (C -m ⊗ m)] = 1 N -1 Ĵ2 -j(j + 1)N - 1 N Ĵ 2 . ( 57
)
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:

Tr[N P (C -m ⊗ m)] = 1 N -1 ( Ĵz ) 2 -Qz - 1 N Ĵz 2 . ( 58
)
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. [START_REF] Frérot | Quantum critical metrology[END_REF], 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), na,+1 (i), na,-1 (i), na,-2 (i)} a∈{x,y,z} . This choice leads to the Tr(C) term:

Tr(C) = 1 N (N -1) a∈{x,y,z} s =0 [ Na,s ] 2 -Na,s (59)
1 N (N -1) a∈{x,y,z}   s =0 [ Na,s ] 2 -(N -Na,0   (60) = 1 N (N -1) a∈{x,y,z}   s =0 [ Na,s ] 2 + Na,0   - 3 N -1 (61)
The method reveals the optimal projector P x = P y

= P z = 1 2       1 0 0 -1 0 1 -1 0 0 -1 1 0 -1 0 0 1      
. This yields the terms:

Tr[N P C] = 1 2(N -1)   a∈{x,y,z} [ Na,+2 -Na,-2 ] 2 + [ Na,+1 -Na,-1 ] 2 + Na,0 -3N   (62) Tr[N P (m ⊗ m)] = 1 2N a∈{x,y,z} Na,+2 -Na,-2 2 + Na,+1 -Na,-1 2 . ( 63
)
We introduce the notation σa,s = Na,s -Na,-s , the population imbalance between opposed levels along direction a. Together with the notation 1 = (1, 1, 1), we can rewrite these equations as:

2 + [ σ2 ] 2 - N -1 N [ σ1 2 + σ2 2 ] + 1 - 2 N N0 • 1 - 2 N s =0 [ Ns ] 2 + 3 , (67)
namely, Eq. [START_REF] Frérot | Quantum critical metrology[END_REF].

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 -na,-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), na,-1 (i)} K a={x,y,z} } N i=1 , where na,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), na,0 (i)} K a={x,y,z} } N i=1 .

Figure 1 :

 1 Figure 1: (a) Comparison of the population-based witness [Eq. (6)] (solid line) with the SSI Var( Ĵ) ≥ N j (dashed line) for spin-1 states (N = 100 atoms). The light grey area are data points which cannot be achieved (see text and Appendix E). The dotted line indicates the ground state of Eq. (5) for both antiferromagnetic interactions (c = 1), varying the quadratic Zeeman term (-10 ≤ q ≤ 10). Three representative states are also indicated: the polar state |0, N, 0 (q → ∞) which is separable, the spin singlet s.t. Ĵ2 |Ψ = 0 (q = 0), and the twin-Fock state |N/2, 0, N/2 (q → -∞). (b) Comparison of the population-based witness [Eq.[START_REF] Duan | Quantum entanglement in spinor boseeinstein condensates[END_REF]] with a related SSI (see text) for spin-1 states. Graphical conventions as in panel (a). Here, the ground state is computed for ferromagnetic interactions (c = -1). At q = 0, we obtain the Dicke state Ĵ2 |Dicke = N (N + 1)|Dicke , Ĵz |Ψ = 0.

Figure 2 :

 2 Figure 2: Comparison of the population-based entanglement witness of Eq. (9) (solid line) with the related SSI Var( Ĵ) ≥ N j (dashed line) for spin-2 ensembles (N = 24). The light grey area are data points which cannot be achieved (see text and Appendix E). Data in the ground state of Eq. (8) are computed for c = 1, and p = 1 (dotted line) or p = -1 (dashed-dotted line). Four representative states are indicated: the polar state |0, 0, N, 0, 0 (with J 2 /(2N ) = 3), the twin-Fock state |N/2, 0, 0, 0, N/2 and two forms of spin singlets (Var(J) = 0), obtained either as a pair condensate or a triplet condensate (see text).

2 - 1 γ=1f

 21 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: λα λβλβ λα = i d αβγ λγ , (48) are characterized by the structure constants f αβγ . The structure constants are the same for the generators of the SU (d) group with the Lie bracket [ λα , λβ ] := λα λβλβ λα . The representations of SU (d) are labelled by partitions of N in d parts as the permutation group of N d-level atoms. In the totally symmetric sector with respect to the exchange of atoms, the matrix representation of the generators is computed { Λα = I[ λα ]} d 2 -1

2(N - 1 )

 1 Tr(C) = 2 N s =0 [ Ns ] 2 + 2 N N0 • 1 -6 (64) 2N (N -1)Tr[P C] = [ σ1 ] 2 + [ σ2 ] 2 + N0 • 1 -Together with the bound β = 3/2 found numerically, we find:2(N -1) [N Tr[P (Cm ⊗ m)] -Tr(C) + β] = [ σ1 ]

  for {n \ n0 } 1.25 1.5 1.5 1.33691.2926

Table 2 :

 2 Numerical separable bound β for the basis of local measurements {{n a,s } j s=-j =0 } a∈{x,y,z} for different spin-j.

Accepted in Quantum 2022-12-21, click title to verify. Published under CC-BY 4.0.

Acknowledgments

We are grateful to Fabrice Gerbier for useful discussions and comments. GM and ML acknowledges support from: ERC AdG NO-QIA; Agencia Estatal de Investigación (R & D project CEX2019-000910-S, funded by MCIN/ AEI/10.13039/501100011033, Plan National FIDEUA PID2019-106901GB-I00, FPI, QUAN-TERA MAQS PCI2019-111828-2, Proyectos de I+D+I "Retos Colaboración" QUSPIN RTC2019-007196-7); Fundació Cellex; Fun-dació Mir-Puig; Generalitat de Catalunya through the European Social Fund FEDER and CERCA program (AGAUR Grant No. 2017 SGR 134, QuantumCAT U16-011424, co-funded by ERDF Operational Program of Catalonia 2014-2020); EU Horizon 2020 FET-OPEN OPTOlogic (Grant No 899794); National Science Centre, Poland (Symfonia Grant No. 2016/20/W/ST4/00314); European Union's Horizon 2020 research and innovation programme under the Marie-Skłodowska-Curie grant agreement No 101029393 (STREDCH) and No 847648 ("La Caixa" Junior Leaders fellowships ID100010434: LCF/BQ/PI19/11690013, LCF/BQ/PI20/11760031, LCF/BQ/PR20/11770012, LCF/BQ/PR21/11840013). IF acknowledges support from the Agence Nationale de la Recherche (Qu-DICE project ANR-PRC-CES47), the John Templeton Foundation (Grant No. 61835), European Union's Horizon 2020 research and innovation programme under the Marie-Skłodowska-Curie grant agreement No 101031549 (QuoMoDys).

Code availability. The code used to diagonalize the spin-1 and spin-2 Hamiltonians, to compute the collective populations in the ground state, and to infer the optimal entanglement witnesses presented in this paper, is available at https://github.com/GuillemMRR/EW_ Zeeman. See Appendix H for a concrete example.

of:

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 of L, 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), [START_REF] Duan | Quantum entanglement in spinor boseeinstein condensates[END_REF] and their generalization to arbitrary integer spins Eqs. (9), [START_REF] Evrard | Observation of fragmentation of a spinor bose-einstein condensate[END_REF] 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 projectors na,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. (3). We evaluated this bound numerically, and the corresponding values are: Table 1: Numerical separable bound β for the basis of local measurements {{ Ŝa (i), na,0 (i)} a={x,y,z} } i={1,2,..,N } for different spin-j.

In the specific case of j = 1, the EW inequalities can be simplified using the identities: a∈{x,y,z} ( Na,+1 + Na,-1 ) = 3N -a∈{x,y,z} Na,0 = N i=1 [ Ŝ(i)] 2 = j(j + 1)N = 2N . Together with the bound β j=1 = 3/2, we obtain Eqs. (6), [START_REF] Duan | Quantum entanglement in spinor boseeinstein condensates[END_REF] from respectively Eqs. (52) and (53).

Derivation of Tr[C].

To derive these EWs from our central inequality, Eq. (4), one has to compute Tr[C], which reads:

where Ĵ = ( Ĵx , Ĵy , Ĵz ) is the collective spin, N0 = ( Nx,0 , Ny,0 , Nz,0 ) and 1 = (1, 1, 1). We used that a [ Ŝa (i)] 2 = j(j +1), and that [n a,0 (i)] 2 = na,0 (i). The two EW are obtained by choosing two different projectors in Eq. (4). Recall that, since correlations in different orientations are not observed, P has a block-diagonal form (P x , P y , P z ).

H.2 Inferring the necessary statistics from the experiment

According to the discussion of the main text, we need to compute the following expectation values:

Na,+1 ( Na,+1 -1) Na,+1 Na,-1 Na,-1 Na,+1 Na,-1 ( Na,-1 -1)

a∈{x,y,z} , (68) where Na,s = N i=1 na,s (i), that is, the population on hyperfine level s along direction a. In our numerical illustration, we can consider N = 100 atoms in the ground state of the Hamiltonian of Eq. ( 5) for c = q = 1. We compute the expectation values of Eq. (68) following the procedure detailed in Appendix D, which yields: m x = m y = (0.499, 0.499), m z = (0.002, 0.002) , C x = C y = 0.248 0.251 0.251 0.248 , C z = 0.000 0.000 0.000 0.000 .

(69) One can prepare such state on an unpolarized spin-1 BEC adiabatically by ramping down the magnetic field to ∼ √ q from q/c 1. Then, one needs to perform a collective Stern-Gerlach splitting along three orthogonal spatial directions, as achieved e.g. in ref. [START_REF] Evrard | Observation of fragmentation of a spinor bose-einstein condensate[END_REF]. The statistics of (68) can readily be inferred by counting the number of atoms on each magnetic sublevel via fluorescence imaging [START_REF] Evrard | Observation of fragmentation of a spinor bose-einstein condensate[END_REF].

H.3 Detecting quantum entanglement

Now we obtain the matrices {C a -m a ⊗ m a := Ca } a∈{x,y,z} , which are central to our method. With the values in Eq. (69), we have: 2 in Appendix G), we have W + β = -0.300 < 0. Therefore, according to the main inequality Eq. (4) entanglement is detected. That is, there is no separable state compatible with the statistics contained in {C a , m a } a∈{x,y,,z} .

For the present case Eq. (69), the projectors to the negative spectrum are:

Remarkably, we find the same optimal projectors regardless the value of the quadratic Zeeman strength q. Finally, from the set of local measurements, the projectors {P x , P y , P z } and β one can derive the entanglement for all N in terms of collective spin and population to the zero Zeeman sublevel as done e.g. in Appendix F. Once the entanglement witness Eq. ( 6) is derived, its violation can be verified equivalently from the expectation values Ĵ2 , Ĵ , N2 0 .