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

Guillem Müller-Rigat<EMAIL_ADDRESS>ICFO - Institut de Ciencies
Fotoniques, The Barcelona Institute of Science and Technology, 08860
Castelldefels (Barcelona), Spain  Maciej Lewenstein  ICFO - Institut de
Ciencies Fotoniques, The Barcelona Institute of Science and Technology, 08860
Castelldefels (Barcelona), Spain  ICREA, Pg. Lluís Companys 23, 08010
Barcelona, Spain  Irénée Frérot<EMAIL_ADDRESS>Univ Grenoble
Alpes, CNRS, Grenoble INP, Institut Néel, 38000 Grenoble, France

###### Abstract

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,\dots,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 $d$-level quantum systems.

##  1 Introduction

To prepare and detect quantum-entangled states has become a central goal for
experimental many-body 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 [31, 17], 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 so-called
spin-squeezing inequalities (SSI), which involve measurements of first- and
second-moments of collective spin operators [22, 44, 37, 36, 39, 41, 29, 42,
31]. SSI and their generalizations [39, 41, 42, 16] have found an impressively
broad range of successful applications, from detecting entanglement in cold
and hot atomic ensembles [31, 23] to unveiling many-body Bell non-locality
[35, 8, 32, 30], from probing quantum phase transitions [11] and quantum
quenches [6, 5] 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 [13, 33, 9, 10, 18, 4,
43, 19, 28, 45, 24, 27, 26, 25, 40, 3, 1, 14], mixtures [20, 2], or effective
qudits ensembles [34]. In these systems, the populations of all Zeeman
sublevels $s=-j,-j+1,\dots,j$ can be measured, typically by fluorescence
imaging after collective spin rotations and Stern-Gerlach splitting [18, 45,
33, 9, 10, 28, 24, 25, 26, 14]. 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 [18, 45,
33, 9, 10]. 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 $\hat{\rho}$, the many-
body density operator of the system, as a statistical mixture of product
states over individual atoms: $\hat{\rho}\neq\hat{\rho}_{\rm
sep}:=\sum_{\lambda}p_{\lambda}\otimes_{i=1}^{N}\hat{\rho}_{\lambda}(i)$.
Here, $\hat{\rho}_{\lambda}(i)$ is an arbitrary internal state (pure or mixed)
for atom $i=1,\dots,N$, and $p_{\lambda}>0$ is an arbitrary probability
distribution, defining a separable state $\hat{\rho}_{\rm sep}$. Our starting
point is to consider for each atom a set of $K$ different observables
$\\{\hat{o}_{k}(i)\\}_{k=1}^{K}$, some of which must be non-commuting in order
to potentially detect entanglement. We then assume that data of the following
form are available:

$\displaystyle
m_{k}=\frac{1}{N}\sum_{i=1}^{N}\langle\hat{o}_{k}(i)\rangle\leavevmode\nobreak\
,$ (1) $\displaystyle C_{kl}=\frac{1}{N(N-1)}\sum_{i\neq
j}\langle\hat{o}_{k}(i)\hat{o}_{l}(j)\rangle\leavevmode\nobreak\ ,$ (2)

where $\langle\dots\rangle={\rm Tr}(\hat{\rho}\dots)$ 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 $\beta$ defined as:

$\beta=\frac{1}{N}\sum_{i=1}^{N}\max_{|\psi\rangle}\left\\{\sum_{k=1}^{K}|\langle\psi|\hat{o}_{k}(i)|\psi\rangle|^{2}\right\\}\leavevmode\nobreak\
,$ (3)

where $|\psi\rangle$ is an arbitrary internal state of atom $i$. Depending on
the choice of observables $\hat{o}_{k}(i)$, the bound $\beta$ 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 $\hat{o}_{k}(i)$ are the same for all atoms (they are
the components of the spin, or projectors onto spin states), and we simply
parametrize $|\psi\rangle=\sum_{\alpha=1}^{d}c_{\alpha}|\alpha\rangle$ with
its complex coefficients $c_{\alpha}$ in a fixed basis.
$\sum_{k=1}^{K}|\langle\psi|\hat{o}_{k}(i)|\psi\rangle|^{2}$ is then a
function of the coefficient $c_{\alpha}$ (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\leavevmode\nobreak\ {\rm Tr}[P(C-{\bf m}\otimes{\bf m})]\geq{\rm
Tr}(C)-\beta\leavevmode\nobreak\ ,$ (4)

where we introduced the notation $[{\bf m}\otimes{\bf m}]_{kl}=m_{k}m_{l}$,
and $P$ is any projector (namely a $K\times 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-{\bf m}\otimes{\bf
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 populations $\hat{N}_{a,s}$
of all Zeeman sublevels $s=-j,\dots,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 $\langle\hat{N}_{a,s}\rangle$ and their correlations
$\langle\hat{N}_{a,s}\hat{N}_{a,s^{\prime}}\rangle$ can be inferred in this
manner. In the general framework presented above, the single-atom observables
correspond to the projectors onto the spin states
$\hat{n}_{a,s}(i)=|s_{a}\rangle_{i}\langle 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}=\langle\hat{N}_{a,s}\rangle/N$, and the diagonal
blocks of the correlation matrix corresponding to a fixed orientation
$C_{(a,s)(a,s^{\prime})}=\langle\hat{N}_{a,s}[\hat{N}_{a,s^{\prime}}-\delta_{s,s^{\prime}}]\rangle/[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^{\prime})}-m_{a,s}m_{a,s^{\prime}}$ for the different
quantization axes $a$.

##  2 Applications to a spin-1 BEC

We first illustrate our approach on a BEC of spin-1 atoms, as investigated
recently in several experiments [18, 45, 33, 9, 10]. 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\rangle=\otimes_{i=1}^{N}|0_{z}\rangle_{i}$ (with
$|N_{-},N_{0},N_{+}\rangle$ 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):

$\hat{H}_{j=1}=\frac{c}{N}\hat{\mathbf{J}}^{2}+q\hat{Q}_{z}\leavevmode\nobreak\
,$ (5)

where $\hat{J}_{a}=\sum_{s=-j}^{j}s\hat{N}_{a,s}$ is the collective spin in
direction $a\in\\{x,y,z\\}$ and
$\hat{\mathbf{J}}^{2}=\hat{J}_{x}^{2}+\hat{J}_{y}^{2}+\hat{J}_{z}^{2}$;
$\hat{Q}_{z}=\sum_{s=-j}^{j}s^{2}\hat{N}_{z,s}$ is the collective quadrupole
operator; and $\\{x,y,z\\}$ is an orthornormal basis of $\mathbb{R}^{3}$.
Notice that the linear Zeeman term, proportional to $\hat{J}_{z}$, is omitted
as it commutes with $\hat{\mathbf{J}}^{2}$ and acts trivially on the initial
state: $\hat{J}_{z}|0,N,0\rangle=0$. The initial state is also the ground
state of $\hat{Q}_{z}$; and by varying the magnetic field intensity ($q$) one
can adiabatically prepare the ground state of $\hat{H}_{j=1}$ for either
ferromagnetic interactions ($c\leq 0$) or antiferromagnetic interactions
($c\geq 0$). Notice that in the case of ferromagnetic interactions, two
quantum phase transitions occur at $q/|c|=\pm 4$, so that adiabaticity is
possible only for finite $N$. As the Zeeman sublevel populations are related
by $\sum_{s=-j}^{j}\hat{N}_{a,}=N$ for all directions, $m_{a,s}$ and
$C_{(a,s)(a,s^{\prime})}$ for all $s,s^{\prime}$ are not all independent. As
local observables in Eqs. (1) and (2) we used the spin
$\hat{S}_{a}(i)=\sum_{s=-j}^{j}s\hat{n}_{a,s}(i)$ and the projector onto the
Zeeman sublevel $s=0$, $\hat{n}_{a,0}(i)$. We use these data as input to Eq.
(4), with bound $\beta_{j=1}=3/2$ in Eq. (3) (see Appendix F). Notice that for
$j=1$, we have $\hat{Q}_{a}=\hat{N}_{a,+}+\hat{N}_{a,-}=N-\hat{N}_{a,0}$.

[width=]spin_1.png

Figure 1: (a) Comparison of the population-based witness [Eq. (6)] (solid
line) with the SSI ${\rm Var}(\hat{\bf J})\geq Nj$ (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\leq q\leq 10$). Three representative states are
also indicated: the polar state $|0,N,0\rangle$ ($q\to\infty$) which is
separable, the spin singlet s.t. $\hat{\bf J}^{2}|\Psi\rangle=0$ ($q=0$), and
the twin-Fock state $|N/2,0,N/2\rangle$ ($q\to-\infty$). (b) Comparison of the
population-based witness [Eq. (7)] 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 $\hat{\bf J}^{2}|{\rm Dicke}\rangle=N(N+1)|{\rm Dicke}\rangle$,
$\hat{J}_{z}|\Psi\rangle=0$.

#### Antiferromagnetic interactions.

For antiferromagnetic interactions ($c>0$) and $q/c\lesssim 4.4$ ($N=100$), we
find violation of Eq. (4), with an entanglement witness (EW) which reads (see
Appendix F):

$\mathrm{Var}(\hat{\mathbf{J}})-\frac{1}{N-1}\langle\hat{\mathbf{N}}_{0}^{2}\rangle\geq\frac{N(N-3)}{2(N-1)}\
,$ (6)

with the notations
$\mathrm{Var}(\hat{\mathbf{J}})=\mathrm{Var}(\hat{J}_{x})+\mathrm{Var}(\hat{J}_{y})+\mathrm{Var}(\hat{J}_{z})$,
$\mathrm{Var}(\hat{O})=\langle\hat{O}^{2}\rangle-\langle\hat{O}\rangle^{2}$
and $\hat{\mathbf{N}}_{0}=(\hat{N}_{x,0},\hat{N}_{y,0},\hat{N}_{z,0})$. The EW
of Eq. (6) is maximally violated at $q=0$, where the ground state of Eq. (5)
is a spin singlet ($\mathrm{Var}(\hat{\mathbf{J}})=0$). In particular, we also
applied our method using as input the experimental data of ref. [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
$\mathrm{Var}(\hat{\mathbf{J}})\geq Nj$ (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
$\langle\hat{\mathbf{J}}\rangle=\mathbf{0}$ (vanishing mean spin, as is the
case in the ground state of $\hat{H}_{j=1}$, Eq. (5)) in the plane
$(x=\langle\hat{\mathbf{J}}^{2}\rangle/N,y=\langle\hat{\mathbf{N}}_{0}^{2}\rangle/[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
$\hat{H}_{j=1}$ for $c>0$ and both $q>0$ and $q<0$. More generally, not all
combinations of $\langle\hat{\bf J}^{2}\rangle$ and $\langle\hat{\bf
N}_{0}^{2}\rangle$ 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 $\lambda_{1}\hat{\bf
J}^{2}+\lambda_{2}\hat{\bf N}_{0}^{2}$ for all values of
$\lambda_{1},\lambda_{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 $\hat{H}_{j=1}$, three
limiting cases are especially illustrative. For $q\to+\infty$, the ground
state is the polar state $|0,N,0\rangle=\otimes_{i=1}^{N}|0_{z}\rangle_{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 $\mathrm{Var}(\hat{\mathbf{J}})\geq Nj$. In-between, we notice that our EW
detects entanglement in the ground state of $\hat{H}_{j=1}$ for a larger range
of values of $q$ than SSIs. Finally, for $q\to-\infty$, the ground state is a
twin Fock state $|N/2,0,N/2\rangle$. 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=\pm 1$). Yet, the
twin-Fock states robustly violates the EW of Eq. (6) (with
$\mathrm{Var}(\hat{\mathbf{J}})=N$ and
$\langle\hat{\mathbf{N}}_{0}^{2}\rangle=N(3N+2)/4$, so that the l.h.s is $\sim
N/4$, while the r.h.s is $\sim 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|\lesssim 1.2$ ($N=100$), we
find violation of (see Appendix F):

$\displaystyle(N-1){\rm
Var}(\hat{J}_{z})-N\langle\hat{Q}_{z}\rangle-\langle\hat{J}_{x}^{2}+\hat{J}_{y}^{2}\rangle\rangle$
$\displaystyle+\frac{3N(N+1)}{2}\geq\langle\hat{\mathbf{N}}_{0}^{2}\rangle\ .$
(7)

Maximal violation is found close to $q=0$ where the state is a Dicke state,
such that $\hat{\bf J}^{2}|{\rm Dicke}\rangle=Nj(Nj+1)|{\rm Dicke}\rangle$
(maximal total spin), and $\hat{J}_{z}|{\rm Dicke}\rangle=0$. Eq. (7) can be
compared to another SSI of ref. [42] violated by the Dicke state, namely
$(N-1){\rm
Var}(\hat{J}_{z})-N\langle\hat{Q}_{z}\rangle-\langle\hat{J}_{x}^{2}+\hat{J}_{y}^{2}\rangle+Nj(Nj+1)\geq
0$. In Fig. 1(b), we compare their power in detecting entanglement in the
ground state of $\hat{H}_{j=1}$. In particular, for $q\lesssim|c|$, the ground
state of $\hat{H}_{j=1}$ is very close to the twin-Fock state
$|N/2,0,N/2\rangle$, which robustly violates Eq. (7): the l.h.s is equal to
$N(N+1)/2\sim N^{2}/2$, while on the r.h.s
$\langle\hat{\mathbf{N}}_{0}^{2}\rangle=N(3N+2)/4\sim 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\rangle$, stabilized as the ground
state of $\hat{H}_{j=1}$ in the limit $q\ll|c|$; yet no generalized SSI for
spin-1 ensembles is violated by the twin-Fock state.

[width=]spin_2.pdf

Figure 2: Comparison of the population-based entanglement witness of Eq. (9)
(solid line) with the related SSI ${\rm Var}(\hat{\bf J})\geq Nj$ (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\rangle$ (with $\langle\mathbf{J}^{2}/(2N)\rangle=3$), the twin-
Fock state $|N/2,0,0,0,N/2\rangle$ and two forms of spin singlets (${\rm
Var}({\bf J})=0$), obtained either as a pair condensate or a triplet
condensate (see text).

##  3 Applications to a spin-2 BEC

The Hamiltonian of a BEC of $j=2$ atoms interacting via pairwise isotropic
interactions reads (see Appendix C):

$\hat{H}_{j=2}=\frac{c}{N}\hat{\mathbf{J}}^{2}+\frac{p}{N}\hat{\Theta}_{2}^{\dagger}\hat{\Theta}_{2}+q\hat{Q}_{z}\
.$ (8)

In contrast to the spin-1 case [Eq. (5)], the Hamiltonian contains the
additional pairing term $\hat{\Theta}_{2}^{\dagger}\hat{\Theta}_{2}$ where
$\hat{\Theta}_{2}=\sum_{s=-j}^{j}(-1)^{s}\hat{a}_{s}\hat{a}_{-s}$ annihilates
time-reversal pairs $\pm s$, with $\hat{a}_{s}$ annihilating an atom in the
spin state $|s_{z}\rangle$. Notice that
$[\hat{\mathbf{J}}^{2},\hat{\Theta}_{2}^{\dagger}\hat{\Theta}_{2}]=0$. We
explored the ground state of $\hat{H}_{j=2}$ to find novel violated EW. Using
as local operators $\hat{S}_{a}(i)$ (the spin) and $\hat{n}_{a,0}(i)$ (the
projector onto the $s=0$ spin state) for $a\in\\{x,y,z\\}$, the bound in Eq.
(3) is numerically found as $\beta_{2}\approx 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):

$\mathrm{Var}(\hat{\mathbf{J}})-\frac{1}{N-1}\langle\hat{\mathbf{N}}_{0}(\hat{\mathbf{N}}_{0}-\mathbf{1})\rangle\geq
N[j(j+1)-\beta_{j}]\ ,$ (9)

and:

$\displaystyle(N-1){\rm
Var}(\hat{J}_{z})-N\langle\hat{Q}_{z}\rangle-\langle\hat{J}_{x}^{2}+\hat{J}_{y}^{2}\rangle$
$\displaystyle+N[(N-1)\beta_{j}+j(j+1)]\geq\langle\hat{\mathbf{N}}_{0}(\hat{\mathbf{N}}_{0}-\mathbf{1})\rangle\leavevmode\nobreak\
,$ (10)

where $\mathbf{1}=(1,1,1)$ and $\beta_{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\rangle$ (obtained in the limit
$q\to-\infty$) 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 ($\langle\hat{\bf J}^{2}\rangle=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”
$(\hat{\Theta}^{\dagger}_{2})^{N/2}\ket{\mathrm{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”
$(\hat{\Theta}^{\dagger}_{3})^{N/3}\ket{\mathrm{vac.}}$, where
$\Theta^{\dagger}_{3}=-9[(a^{\dagger}_{-1})^{2}a^{\dagger}_{2}+(a^{\dagger}_{1})^{2}a^{\dagger}_{-2}]+\sqrt{6}a^{\dagger}_{0}[-(a^{\dagger}_{0})^{2}+3a^{\dagger}_{1}a^{\dagger}_{-1}+6a^{\dagger}_{2}a^{\dagger}_{-2}]$
is the singlet-trio creation operator. Both states are entangled and violate
Eq. (9). However, considering as local measurements
$\\{\hat{n}_{a,+2}(i),\hat{n}_{a,+1}(i),\hat{n}_{a,-1}(i),\hat{n}_{a,-2}(i)\\}_{a\in\\{x,y,z\\}}$,
we found the following EW violated only by the pair condensate and not by the
triplet condensate (see Appendix G):

$\displaystyle\langle\hat{\boldsymbol{\sigma}}_{2}^{2}+\hat{\boldsymbol{\sigma}}_{1}^{2}\rangle-\frac{N-1}{N}(\langle\hat{\boldsymbol{\sigma}}_{2}\rangle^{2}+\langle\hat{\boldsymbol{\sigma}}_{1}\rangle^{2})$
$\displaystyle-\frac{2}{N}\langle\hat{\mathbf{N}}_{+2}^{2}+\hat{\mathbf{N}}_{+1}^{2}+\hat{\mathbf{N}}_{-1}^{2}+\hat{\mathbf{N}}_{-2}^{2}\rangle$
$\displaystyle+(1-\frac{2}{N})\langle\hat{\mathbf{N}}_{0}\cdot\mathbf{1}\rangle+3\geq
0\ ,$ (11)

where we denoted by $\hat{\sigma}_{a,s}=\hat{N}_{a,s}-\hat{N}_{a,-s}$ the
population imbalance between opposed levels along direction $a$. Notice that
$\\{\hat{\sigma}_{x,s},\hat{\sigma}_{y,s},\hat{\sigma}_{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.

##  4 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 two-body
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 spatially-
structured systems by using different local observables $\hat{o}_{k}(i)$ for
each atom; a straightforward example is to introduce local phases
$e^{i\phi_{k}(i)}$, leaving the bound $\beta$ 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 $\beta$ 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.

#### 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.

We are grateful to Fabrice Gerbier for useful discussions and comments. GM and
ML acknowledges support from: ERC AdG NOQIA; Agencia Estatal de Investigación
(R$\And$D project CEX2019-000910-S, funded by MCIN/ AEI/10.13039/501100011033,
Plan National FIDEUA PID2019-106901GB-I00, FPI, QUANTERA MAQS
PCI2019-111828-2, Proyectos de I+D+I “Retos Colaboración” QUSPIN
RTC2019-007196-7); Fundació Cellex; Fundació 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).

## References

  * Bataille et al. [2020] P. Bataille, A. Litvinov, I. Manai, J. Huckans, F. Wiotte, A. Kaladjian, O. Gorceix, E. Maréchal, B. Laburthe-Tolra, and M. Robert-de Saint-Vincent.  Adiabatic spin-dependent momentum transfer in an su($n$) degenerate fermi gas.  _Phys. Rev. A_ , 102:013317, Jul 2020.  doi: 10.1103/PhysRevA.102.013317.  URL https://link.aps.org/doi/10.1103/PhysRevA.102.013317.
  * Bhatt et al. [2022] Rohit Prasad Bhatt, Jan Kilinc, Lilo Höcker, and Fred Jendrzejewski.  Stochastic dynamics of a few sodium atoms in presence of a cold potassium cloud.  _Scientific Reports_ , 12(1), Feb 2022.  ISSN 2045-2322.  doi: 10.1038/s41598-022-05778-8.  URL http://dx.doi.org/10.1038/s41598-022-05778-8.
  * Chalopin et al. [2018] Thomas Chalopin, Chayma Bouazza, Alexandre Evrard, Vasiliy Makhalov, Davide Dreon, Jean Dalibard, Leonid A. Sidorenkov, and Sylvain Nascimbene.  Quantum-enhanced sensing using non-classical spin states of a highly magnetic atom.  _Nature Communications_ , 9(1), nov 2018.  doi: 10.1038/s41467-018-07433-1.  URL https://doi.org/10.1038%2Fs41467-018-07433-1.
  * Chomaz et al. [2016] L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino.  Quantum-fluctuation-driven crossover from a dilute bose-einstein condensate to a macrodroplet in a dipolar quantum fluid.  _Phys. Rev. X_ , 6:041039, Nov 2016.  doi: 10.1103/PhysRevX.6.041039.  URL https://link.aps.org/doi/10.1103/PhysRevX.6.041039.
  * Comparin et al. [2022a] Tommaso Comparin, Fabio Mezzacapo, and Tommaso Roscilde.  Multipartite entangled states in dipolar quantum simulators.  _Physical Review Letters_ , 129(15), oct 2022a.  doi: 10.1103/physrevlett.129.150503.  URL https://doi.org/10.1103%2Fphysrevlett.129.150503.
  * Comparin et al. [2022b] Tommaso Comparin, Fabio Mezzacapo, and Tommaso Roscilde.  Robust spin squeezing from the tower of states of u(1)-symmetric spin hamiltonians.  _Phys. Rev. A_ , 105:022625, Feb 2022b.  doi: 10.1103/PhysRevA.105.022625.  URL https://link.aps.org/doi/10.1103/PhysRevA.105.022625.
  * Duan et al. [2002] L.-M. Duan, J. I. Cirac, and P. Zoller.  Quantum entanglement in spinor bose-einstein condensates.  _Phys. Rev. A_ , 65:033619, Feb 2002.  doi: 10.1103/PhysRevA.65.033619.  URL https://link.aps.org/doi/10.1103/PhysRevA.65.033619.
  * Engelsen et al. [2017] Nils J. Engelsen, Rajiv Krishnakumar, Onur Hosten, and Mark A. Kasevich.  Bell Correlations in Spin-Squeezed States of 500 000 Atoms.  _Phys. Rev. Lett._ , 118(14):140401, April 2017\.  doi: 10.1103/PhysRevLett.118.140401.  URL https://link.aps.org/doi/10.1103/PhysRevLett.118.140401.
  * Evrard et al. [2021a] Bertrand Evrard, An Qu, Jean Dalibard, and Fabrice Gerbier.  From many-body oscillations to thermalization in an isolated spinor gas.  _Phys. Rev. Lett._ , 126:063401, Feb 2021a.  doi: 10.1103/PhysRevLett.126.063401.  URL https://link.aps.org/doi/10.1103/PhysRevLett.126.063401.
  * Evrard et al. [2021b] Bertrand Evrard, An Qu, Jean Dalibard, and Fabrice Gerbier.  Observation of fragmentation of a spinor bose-einstein condensate.  _Science_ , 373(6561):1340–1343, 2021b.  doi: 10.1126/science.abd8206.  URL https://www.science.org/doi/abs/10.1126/science.abd8206.
  * Frérot and Roscilde [2018] Irénée Frérot and Tommaso Roscilde.  Quantum critical metrology.  _Phys. Rev. Lett._ , 121:020402, Jul 2018.  doi: 10.1103/PhysRevLett.121.020402.  URL https://link.aps.org/doi/10.1103/PhysRevLett.121.020402.
  * Frérot et al. [2021] Irénée Frérot, Flavio Baccari, and Antonio Acín.  Unveiling quantum entanglement in many-body systems from partial information.  _arXiv preprint arXiv:2107.03944_ , 2021.  doi: https://doi.org/10.48550/arXiv.2107.03944.  URL https://arxiv.org/abs/2107.03944.
  * Frisch et al. [2014] Albert Frisch, Michael Mark, Kiyotaka Aikawa, Francesca Ferlaino, John L. Bohn, Constantinos Makrides, Alexander Petrov, and Svetlana Kotochigova.  Quantum chaos in ultracold collisions of gas-phase erbium atoms.  _Nature_ , 507(7493):475–479, Mar 2014.  ISSN 1476-4687.  doi: 10.1038/nature13137.  URL http://dx.doi.org/10.1038/nature13137.
  * Gabardos et al. [2020] Lucas Gabardos, Bihui Zhu, Steven Lepoutre, Ana Maria Rey, Bruno Laburthe-Tolra, and Laurent Vernac.  Relaxation of the collective magnetization of a dense 3d array of interacting dipolar $s=3$ atoms.  _Phys. Rev. Lett._ , 125:143401, Sep 2020.  doi: 10.1103/PhysRevLett.125.143401.  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.143401.
  * Georgescu et al. [2014] I. M. Georgescu, S. Ashhab, and Franco Nori.  Quantum simulation.  _Rev. Mod. Phys._ , 86:153–185, Mar 2014.  doi: 10.1103/RevModPhys.86.153.  URL https://link.aps.org/doi/10.1103/RevModPhys.86.153.
  * Gessner et al. [2019] Manuel Gessner, Augusto Smerzi, and Luca Pezzè.  Metrological nonlinear squeezing parameter.  _Phys. Rev. Lett._ , 122:090503, Mar 2019.  doi: 10.1103/PhysRevLett.122.090503.  URL https://link.aps.org/doi/10.1103/PhysRevLett.122.090503.
  * Gilmore et al. [2021] Kevin A. Gilmore, Matthew Affolter, Robert J. Lewis-Swan, Diego Barberena, Elena Jordan, Ana Maria Rey, and John J. Bollinger.  Quantum-enhanced sensing of displacements and electric fields with two-dimensional trapped-ion crystals.  _Science_ , 373(6555):673–678, Aug 2021.  ISSN 1095-9203.  doi: 10.1126/science.abi5226.  URL http://dx.doi.org/10.1126/science.abi5226.
  * Hamley et al. [2012] C. D. Hamley, C. S. Gerving, T. M. Hoang, E. M. Bookjans, and M. S. Chapman.  Spin-nematic squeezed vacuum in a quantum gas.  _Nature Physics_ , 8(4):305–308, Feb 2012.  ISSN 1745-2481.  doi: 10.1038/nphys2245.  URL http://dx.doi.org/10.1038/nphys2245.
  * Kadau et al. [2016] Holger Kadau, Matthias Schmitt, Matthias Wenzel, Clarissa Wink, Thomas Maier, Igor Ferrier-Barbut, and Tilman Pfau.  Observing the rosensweig instability of a quantum ferrofluid.  _Nature_ , 530(7589):194–197, Feb 2016.  ISSN 1476-4687.  doi: 10.1038/nature16485.  URL http://dx.doi.org/10.1038/nature16485.
  * Kasper et al. [2021] Valentin Kasper, Daniel González-Cuadra, Apoorva Hegde, Andy Xia, Alexandre Dauphin, Felix Huber, Eberhard Tiemann, Maciej Lewenstein, Fred Jendrzejewski, and Philipp Hauke.  Universal quantum computation and quantum error correction with ultracold atomic mixtures.  _Quantum Science and Technology_ , 7(1):015008, 2021.  doi: 10.1088/2058-9565/ac2d39.
  * Kawaguchi and Ueda [2012] Yuki Kawaguchi and Masahito Ueda.  Spinor bose–einstein condensates.  _Physics Reports_ , 520(5):253–381, 2012.  ISSN 0370-1573.  doi: https://doi.org/10.1016/j.physrep.2012.07.005.  URL https://www.sciencedirect.com/science/article/pii/S0370157312002098.  Spinor Bose–Einstein condensates.
  * Kitagawa and Ueda [1993] Masahiro Kitagawa and Masahito Ueda.  Squeezed spin states.  _Phys. Rev. A_ , 47:5138–5143, Jun 1993.  doi: 10.1103/PhysRevA.47.5138.  URL https://link.aps.org/doi/10.1103/PhysRevA.47.5138.
  * Kong et al. [2020] Jia Kong, Ricardo Jiménez-Martínez, Charikleia Troullinou, Vito Giovanni Lucivero, Géza Tóth, and Morgan W. Mitchell.  Measurement-induced, spatially-extended entanglement in a hot, strongly-interacting atomic system.  _Nature Communications_ , 11(1), May 2020.  ISSN 2041-1723.  doi: 10.1038/s41467-020-15899-1.  URL http://dx.doi.org/10.1038/s41467-020-15899-1.
  * Kunkel et al. [2018] Philipp Kunkel, Maximilian Prüfer, Helmut Strobel, Daniel Linnemann, Anika Frölian, Thomas Gasenzer, Martin Gärttner, and Markus K. Oberthaler.  Spatially distributed multipartite entanglement enables epr steering of atomic clouds.  _Science_ , 360(6387):413–416, 2018.  doi: 10.1126/science.aao2254.  URL https://www.science.org/doi/abs/10.1126/science.aao2254.
  * Kunkel et al. [2019] Philipp Kunkel, Maximilian Prüfer, Stefan Lannig, Rodrigo Rosa-Medina, Alexis Bonnin, Martin Gärttner, Helmut Strobel, and Markus K. Oberthaler.  Simultaneous readout of noncommuting collective spin observables beyond the standard quantum limit.  _Phys. Rev. Lett._ , 123:063603, Aug 2019.  doi: 10.1103/PhysRevLett.123.063603.  URL https://link.aps.org/doi/10.1103/PhysRevLett.123.063603.
  * Lange et al. [2018] Karsten Lange, Jan Peise, Bernd Lücke, Ilka Kruse, Giuseppe Vitagliano, Iagoba Apellaniz, Matthias Kleinmann, Géza Tóth, and Carsten Klempt.  Entanglement between two spatially separated atomic modes.  _Science_ , 360(6387):416–418, 2018.  doi: 10.1126/science.aao2035.  URL https://www.science.org/doi/abs/10.1126/science.aao2035.
  * Lucioni et al. [2018] E. Lucioni, L. Tanzi, A. Fregosi, J. Catani, S. Gozzini, M. Inguscio, A. Fioretti, C. Gabbanini, and G. Modugno.  Dysprosium dipolar bose-einstein condensate with broad feshbach resonances.  _Phys. Rev. A_ , 97:060701, Jun 2018.  doi: 10.1103/PhysRevA.97.060701.  URL https://link.aps.org/doi/10.1103/PhysRevA.97.060701.
  * Luo et al. [2017] Xin-Yu Luo, Yi-Quan Zou, Ling-Na Wu, Qi Liu, Ming-Fei Han, Meng Khoon Tey, and Li You.  Deterministic entanglement generation from driving through quantum phase transitions.  _Science_ , 355(6325):620–623, 2017.  doi: 10.1126/science.aag1106.  URL https://www.science.org/doi/abs/10.1126/science.aag1106.
  * Ma et al. [2011] Jian Ma, Xiaoguang Wang, C.P. Sun, and Franco Nori.  Quantum spin squeezing.  _Physics Reports_ , 509(2-3):89–165, dec 2011\.  doi: 10.1016/j.physrep.2011.08.003.  URL https://doi.org/10.1016%2Fj.physrep.2011.08.003.
  * Müller-Rigat et al. [2021] Guillem Müller-Rigat, Albert Aloy, Maciej Lewenstein, and Irénée Frérot.  Inferring nonlinear many-body bell inequalities from average two-body correlations: Systematic approach for arbitrary spin-$j$ ensembles.  _PRX Quantum_ , 2:030329, Aug 2021.  doi: 10.1103/PRXQuantum.2.030329.  URL https://link.aps.org/doi/10.1103/PRXQuantum.2.030329.
  * Pezzè et al. [2018] Luca Pezzè, Augusto Smerzi, Markus K. Oberthaler, Roman Schmied, and Philipp Treutlein.  Quantum metrology with nonclassical states of atomic ensembles.  _Rev. Mod. Phys._ , 90:035005, Sep 2018.  doi: 10.1103/RevModPhys.90.035005.  URL https://link.aps.org/doi/10.1103/RevModPhys.90.035005.
  * Piga et al. [2019] Angelo Piga, Albert Aloy, Maciej Lewenstein, and Irénée Frérot.  Bell correlations at ising quantum critical points.  _Phys. Rev. Lett._ , 123:170604, Oct 2019.  doi: 10.1103/PhysRevLett.123.170604.  URL https://link.aps.org/doi/10.1103/PhysRevLett.123.170604.
  * Qu et al. [2020] An Qu, Bertrand Evrard, Jean Dalibard, and Fabrice Gerbier.  Probing spin correlations in a bose-einstein condensate near the single-atom level.  _Phys. Rev. Lett._ , 125:033401, Jul 2020.  doi: 10.1103/PhysRevLett.125.033401.  URL https://link.aps.org/doi/10.1103/PhysRevLett.125.033401.
  * Ringbauer et al. [2022] Martin Ringbauer, Michael Meth, Lukas Postler, Roman Stricker, Rainer Blatt, Philipp Schindler, and Thomas Monz.  A universal qudit quantum processor with trapped ions.  _Nature Physics_ , 18(9):1053–1057, jul 2022\.  doi: 10.1038/s41567-022-01658-0.  URL https://www.nature.com/articles/s41567-022-01658-0.
  * Schmied et al. [2016] Roman Schmied, Jean-Daniel Bancal, Baptiste Allard, Matteo Fadel, Valerio Scarani, Philipp Treutlein, and Nicolas Sangouard.  Bell correlations in a Bose-Einstein condensate.  _Science_ , 352(6284):441–444, April 2016.  ISSN 0036-8075, 1095-9203.  doi: 10.1126/science.aad8665.  URL http://science.sciencemag.org/content/352/6284/441.
  * Sørensen et al. [2001] A. Sørensen, L.-M. Duan, J. I. Cirac, and P. Zoller.  Many-particle entanglement with bose–einstein condensates.  _Nature_ , 409(6816):63–66, jan 2001.  doi: 10.1038/35051038.  URL https://doi.org/10.1038%2F35051038.
  * Sørensen and Mølmer [2001a] Anders S. Sørensen and Klaus Mølmer.  Entanglement and extreme spin squeezing.  _Phys. Rev. Lett._ , 86:4431–4434, May 2001a.  doi: 10.1103/PhysRevLett.86.4431.  URL https://link.aps.org/doi/10.1103/PhysRevLett.86.4431.
  * Sørensen and Mølmer [2001b] Anders S. Sørensen and Klaus Mølmer.  Entanglement and extreme spin squeezing.  _Phys. Rev. Lett._ , 86:4431–4434, May 2001b.  doi: 10.1103/PhysRevLett.86.4431.  URL https://link.aps.org/doi/10.1103/PhysRevLett.86.4431.
  * Tóth et al. [2009] Géza Tóth, Christian Knapp, Otfried Gühne, and Hans J. Briegel.  Spin squeezing and entanglement.  _Phys. Rev. A_ , 79:042334, Apr 2009.  doi: 10.1103/PhysRevA.79.042334.  URL https://link.aps.org/doi/10.1103/PhysRevA.79.042334.
  * Trautmann et al. [2018] A. Trautmann, P. Ilzhöfer, G. Durastante, C. Politi, M. Sohmen, M. J. Mark, and F. Ferlaino.  Dipolar quantum mixtures of erbium and dysprosium atoms.  _Phys. Rev. Lett._ , 121:213601, Nov 2018.  doi: 10.1103/PhysRevLett.121.213601.  URL https://link.aps.org/doi/10.1103/PhysRevLett.121.213601.
  * Vitagliano et al. [2011] Giuseppe Vitagliano, Philipp Hyllus, Iñigo L. Egusquiza, and Géza Tóth.  Spin squeezing inequalities for arbitrary spin.  _Phys. Rev. Lett._ , 107:240502, Dec 2011.  doi: 10.1103/PhysRevLett.107.240502.  URL https://link.aps.org/doi/10.1103/PhysRevLett.107.240502.
  * Vitagliano et al. [2014] Giuseppe Vitagliano, Iagoba Apellaniz, Iñigo L. Egusquiza, and Géza Tóth.  Spin squeezing and entanglement for an arbitrary spin.  _Phys. Rev. A_ , 89:032307, Mar 2014.  doi: 10.1103/PhysRevA.89.032307.  URL https://link.aps.org/doi/10.1103/PhysRevA.89.032307.
  * Wenzel et al. [2018] Matthias Wenzel, Fabian Böttcher, Jan-Niklas Schmidt, Michael Eisenmann, Tim Langen, Tilman Pfau, and Igor Ferrier-Barbut.  Anisotropic superfluid behavior of a dipolar bose-einstein condensate.  _Physical Review Letters_ , 121(3), Jul 2018.  ISSN 1079-7114.  doi: 10.1103/physrevlett.121.030401.  URL http://dx.doi.org/10.1103/PhysRevLett.121.030401.
  * Wineland et al. [1994] D. J. Wineland, J. J. Bollinger, W. M. Itano, and D. J. Heinzen.  Squeezed atomic states and projection noise in spectroscopy.  _Phys. Rev. A_ , 50:67–88, Jul 1994.  doi: 10.1103/PhysRevA.50.67.  URL https://link.aps.org/doi/10.1103/PhysRevA.50.67.
  * Zou et al. [2018] Yi-Quan Zou, Ling-Na Wu, Qi Liu, Xin-Yu Luo, Shuai-Feng Guo, Jia-Hao Cao, Meng Khoon Tey, and Li You.  Beating the classical precision limit with spin-1 dicke states of more than 10,000 atoms.  _Proceedings of the National Academy of Sciences_ , 115(25):6381–6385, Jun 2018.  ISSN 1091-6490.  doi: 10.1073/pnas.1715105115.  URL http://dx.doi.org/10.1073/pnas.1715105115.

##  Appendix 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:

$\hat{\rho}_{\rm
sep}=\sum_{\lambda}p_{\lambda}\otimes_{i=1}^{N}\hat{\rho}_{\lambda}(i)\leavevmode\nobreak\
,$ (12)

and to arrive at a contradiction. We introduce the notations:

$\displaystyle m_{k}(\lambda,i):={\rm
Tr}[\hat{\rho}_{\lambda}(i)\hat{o}_{k}(i)]$ (13) $\displaystyle
m_{k}(\lambda):=\frac{1}{N}\sum_{i=1}^{N}m_{k}(\lambda,i)\leavevmode\nobreak\
.$ (14)

The key insight of our approach is that, if the state is separable, we may
interpret the quantities $m_{k}(\lambda,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:

$\frac{1}{N}\sum_{i=1}^{N}\sum_{k}[m_{k}(\lambda,i)]^{2}\leq\beta$ (15)

for all $\lambda$. With the above notations, the one-body quantum data [Eq.
(1)] read:

$m_{k}=\frac{1}{N}\sum_{i=1}^{N}\langle\hat{o}_{k}(i)\rangle=\sum_{\lambda}p_{\lambda}m_{k}(\lambda)$
(16)

We then introduce the vector notation ${\bf m}:=(m_{k})$, and the outer
product notation $[{\bf u}\otimes{\bf v}]_{ab}=u_{a}v_{b}$. With these
notations, the two-body quantum data [Eq. (2)] read:

$C=\frac{1}{N(N-1)}\sum_{i\neq j}\sum_{\lambda}p_{\lambda}\leavevmode\nobreak\
{\bf m}(\lambda,i)^{*}\otimes{\bf m}(\lambda,j)\leavevmode\nobreak\ .$ (17)

We now introduce the fluctuations of the $m_{k}(\lambda,i)$ variables:

$\displaystyle A:=\sum_{\lambda}p_{\lambda}[{\bf m}(\lambda)-{\bf
m}]^{*}\otimes[{\bf m}(\lambda)-{\bf m}]\succeq 0$ (18) $\displaystyle
a:=\sum_{\lambda}p_{\lambda}\frac{1}{N}\sum_{i}[{\bf m}(\lambda,i)-{\bf
m}(\lambda)]^{*}\otimes[{\bf m}(\lambda,i)-{\bf m}(\lambda)]\succeq 0$ (19)

The notation $M\succeq 0$ means the the matrix $M$ is positive semidefinite
(PSD), that is, for any vector ${\bf u}$, ${\bf u}^{\dagger}M{\bf u}\geq 0$.
This PSD property is straightforward to verify, since $A$ and $a$ are both of
the form $\sum_{x}w_{x}{\bf v}{\bf v}^{\dagger}$ with $w_{x}\geq 0$ and ${\bf
v}$ a one-dimensional vector.

The matrix $C$ is then:

$\displaystyle C$ $\displaystyle=\frac{N}{N-1}[{\bf m}^{*}\otimes{\bf
m}+A]-\frac{1}{N(N-1)}\sum_{\lambda}p_{\lambda}\sum_{i}{\bf
m}(\lambda,i)^{*}\otimes{\bf m}(\lambda,i)$ (20) $\displaystyle={\bf
m}^{*}\otimes{\bf m}+A-\frac{1}{N-1}a\leavevmode\nobreak\ .$ (21)

From Eq. (20) and Eq. (15)], we have:

${\rm Tr}(C)\geq\frac{N}{N-1}{\rm Tr}[{\bf m}^{*}\otimes{\bf
m}+A]-\frac{\beta}{N-1}\leavevmode\nobreak\ .$ (22)

We then introduce a projector $P$, and write:

$\displaystyle{\rm Tr}[{\bf m}^{*}\otimes{\bf m}+A]$ $\displaystyle={\rm
Tr}[P({\bf m}^{*}\otimes{\bf m}+A)]+{\rm Tr}[(1-P)({\bf m}^{*}\otimes{\bf
m}+A)]$ (23) $\displaystyle={\rm Tr}[P({\bf m}^{*}\otimes{\bf m}+A)]+{\rm
Tr}[(1-P)(C+a/(N-1))]$ (24)

where we used Eq. (21). From the positivity of $A$ and $P$, we have ${\rm
Tr}(PA)\geq 0$. Similarly, from the positivity of $a$ and $1-P$ (which is also
a projector), we have ${\rm Tr}[(1-P)a]\geq 0$. Hence:

${\rm Tr}[{\bf m}^{*}\otimes{\bf m}+A]\geq{\rm Tr}[P({\bf m}^{*}\otimes{\bf
m})]+{\rm Tr}[(1-P)C]$ (25)

Introducing this inequality into Eq. (22), and reorganizing the terms, this
leads to our central result:

$L(P):=N\leavevmode\nobreak\ {\rm Tr}[P(C-{\bf m}^{*}\otimes{\bf m})]\geq{\rm
Tr}(C)-\beta\leavevmode\nobreak\ ,$ (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. [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-{\bf m}^{*}\otimes{\bf m}$.

##  Appendix B Recovering the generalized spin-squeezing inequalities of ref.
[42].

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:
$\hat{J}_{x},\hat{J}_{y},\hat{J}_{z}$. We write here these inequalities (Eq.
(9) of ref. [42]) for completeness:

$\displaystyle\langle\hat{J}_{x}^{2}+\hat{J}_{y}^{2}+\hat{J}_{z}^{2}\rangle$
$\displaystyle\leq$ $\displaystyle Nj(Nj+1)$ (27) $\displaystyle{\rm
Var}(\hat{J}_{x})+{\rm Var}(\hat{J}_{y})+{\rm Var}(\hat{J}_{z})$
$\displaystyle\geq$ $\displaystyle Nj$ (28)
$\displaystyle\langle\hat{J}_{l}^{2}\rangle-\sum_{i=1}^{N}\langle[\hat{S}_{l}(i)]^{2}\rangle+\langle\hat{J}_{m}^{2}\rangle-\sum_{i=1}^{N}\langle[\hat{S}_{m}(i)]^{2}\rangle$
$\displaystyle\leq$ $\displaystyle(N-1)\left[{\rm
Var}(\hat{J}_{k})-\sum_{i=1}^{N}\langle[\hat{S}_{k}(i)]^{2}\rangle\right]$
(29) $\displaystyle+N(N-1)j^{2}$
$\displaystyle\langle\hat{J}_{m}^{2}\rangle-\sum_{i=1}^{N}\langle[\hat{S}_{m}(i)]^{2}\rangle-N(N-1)j^{2}$
$\displaystyle\leq$ $\displaystyle(N-1)\left[{\rm
Var}(\hat{J}_{k})-\sum_{i=1}^{N}\langle[\hat{S}_{k}(i)]^{2}\rangle\right.$
(30) $\displaystyle\left.+{\rm
Var}(\hat{J}_{l})-\sum_{i=1}^{N}\langle[\hat{S}_{l}(i)]^{2}\rangle\right]$

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 $\hat{o}_{k}(i)$ the three spin observables
$\hat{S}_{x}(i),\hat{S}_{y}(i),\hat{S}_{z}(i)$. We then have:

$\displaystyle m_{k}$ $\displaystyle=$
$\displaystyle\langle\hat{J}_{k}\rangle/N$ (31) $\displaystyle C_{kl}$
$\displaystyle=$
$\displaystyle\frac{1}{N(N-1)}\left[\langle\hat{J}_{k}\hat{J}_{l}\rangle-\sum_{i}\langle\hat{S}_{k}(i)\hat{S}_{l}(i)\rangle\right]\leavevmode\nobreak\
.$ (32)

Furthermore, the $\beta$ bound is $\beta=j^{2}$, since for any state of a
spin-$j$ atom, one has $\langle S_{x}(i)\rangle^{2}+\langle
S_{y}(i)\rangle^{2}+\langle S_{z}(i)\rangle^{2}\leq j^{2}$, with equality if
the atomic spin is fully polarized along some direction. Hence we have:

${\rm Tr}(C)-\beta=\frac{1}{N(N-1)}(\langle
J_{x}^{2}+J_{y}^{2}+J_{z}^{2}\rangle)-\frac{1}{N-1}j(j+1)-j^{2}$ (33)

where we used that
$\langle[S_{x}(i)]^{2}+[S_{y}(i)]^{2}+[S_{z}(i)]^{2}\rangle=j(j+1)$. As we
will show explicitly, the 8 SSI [Eqs. (27) to (30)] are obtained by choosing
the projector $P$ in Eq. (26) diagonal, leading to inequalities of the form:

$N\sum_{k\in E}[C_{kk}-m_{k}^{2}]\geq{\rm Tr}(C)-j^{2}$ (34)

with $E$ any subset of $\\{x,y,z\\}$. This last inequality can be written more
explicitly as:

$\sum_{k\in
E}\left[N\langle\hat{J}_{k}^{2}\rangle-(N-1)\langle\hat{J}_{k}\rangle^{2}-N\sum_{i}\langle[S_{k}(i)]^{2}\rangle\right]\geq\langle
J_{x}^{2}+J_{y}^{2}+J_{z}^{2}\rangle-Nj(Nj+1)\leavevmode\nobreak\ .$ (35)

Choosing $E=\emptyset$ directly gives Eq. (27).
Choosing $E=\\{x,y,z\\}$ gives:

$(N-1)\left[{\rm Var}(\hat{J}_{x})+{\rm Var}(\hat{J}_{y})+{\rm
Var}(\hat{J}_{z})\right]\geq Nj(j+1)-Nj(Nj+1)\leavevmode\nobreak\ ,$ (36)

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

$\displaystyle\sum_{k\in
E}\left[N\langle\hat{J}_{k}^{2}\rangle-(N-1)\langle\hat{J}_{k}\rangle^{2}-N\sum_{i}\langle[S_{k}(i)]^{2}\rangle\right]$
$\displaystyle\geq$ $\displaystyle\langle
J_{x}^{2}+J_{y}^{2}+J_{z}^{2}\rangle-\sum_{i}\left\langle[\hat{S}_{x}(i)]^{2}+[\hat{S}_{y}(i)]^{2}+[\hat{S}_{z}(i)]^{2}\right\rangle$
(37) $\displaystyle-N(N-1)j^{2}\leavevmode\nobreak\ .$

Taking then e.g. $E=\\{x\\}$, we find:

$(N-1)\left[{\rm
var}(\hat{J}_{x})-\sum_{i}\langle[\hat{S}_{x}(i)]^{2}\rangle\right]\geq\langle
J_{y}^{2}+J_{z}^{2}\rangle-\sum_{i}\left\langle[\hat{S}_{y}(i)]^{2}+[\hat{S}_{z}(i)]^{2}\right\rangle-N(N-1)j^{2}\leavevmode\nobreak\
,$ (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. [42] are
incorporated in the framework presented in this paper when choosing the spin
components as local observables.

##  Appendix 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 $\mathbb{R}^{3}$. In the single-mode approximation, its
Hamiltonian reads:

$\hat{H}_{0}(\\{g_{f}\\}_{f=0}^{j})=\sum_{s_{1}s_{2},s_{3}s_{4}=-j}^{j}g_{f}C_{s_{1}s_{2},s_{3}s_{4}}^{j\otimes
j=2f}\hat{a}^{\dagger}_{s_{1}}\hat{a}^{\dagger}_{s_{2}}\hat{a}_{s_{3}}\hat{a}_{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_{s_{1}s_{2},s_{3}s_{4}}^{j\otimes
j=2f}=\bra{j,s_{1}}_{j}\otimes\bra{j,s_{2}}_{j}\mathcal{P}^{j\otimes
j=2f}\ket{j,s_{3}}_{j}\otimes\ket{j,s_{4}}_{j}$, with $\mathcal{P}^{j\otimes
j=2f}=\sum_{S=-2f}^{2f}\ket{2f,S}_{j\otimes j}\bra{2f,S}_{j\otimes j}$, where
$S$ is the spin projection of the composite spin-$2f$. The operators
$\hat{a}_{s}(\hat{a}^{\dagger}_{s})$ annihilate (create) a bosonic mode with
spin projection $s\in\\{-j,-j+1,..,j\\}$. The total number of atoms
$\hat{N}=\sum_{s=-j}^{j}\hat{N}_{s}$, where
$\hat{N}_{s}=\hat{a}^{\dagger}_{s}\hat{a}_{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:

$\hat{O}=\sum_{s_{L},s_{R}=-j}^{j}\bra{s_{L}}\hat{o}\ket{s_{R}}\hat{a}^{\dagger}_{s_{L}}\hat{a}_{s_{R}}\
,$ (40)

where $\hat{O}=\sum_{i=1}^{N}{\hat{o}}(i)$ is the collective operator
corresponding to the local observable $\hat{o}$ acting on a single-atom
spin-$j$ degree of freedom. Note that in particular the population operator
$\hat{N}_{s}$ is a collective operator with
$\bra{s_{L}}\hat{n}_{s}\ket{s_{R}}:=\innerproduct{s_{L}}{s}\innerproduct{s}{s_{R}}=\delta_{s_{L},s}\delta_{s,s_{R}}$.

#### Spin 1.

For spin-1, the spin projection onto the cartesian basis $\\{x,y,z\\}$ reads:

$\hat{J}_{x}=\frac{1}{\sqrt{2}}\begin{pmatrix}\hat{a}_{+1}&\hat{a}_{0}&\hat{a}_{-1}\\\
\end{pmatrix}\begin{pmatrix}0&1&0\\\ 1&0&1\\\ 0&1&0\\\
\end{pmatrix}\begin{pmatrix}\hat{a}_{+1}\\\ \hat{a}_{0}\\\
\hat{a}_{-1}\end{pmatrix}$ (41)
$\hat{J}_{y}=\frac{1}{\sqrt{2}}\begin{pmatrix}\hat{a}_{+1}&\hat{a}_{0}&\hat{a}_{-1}\\\
\end{pmatrix}\begin{pmatrix}0&-\mathbb{i}&0\\\ \mathbb{i}&0&-\mathbb{i}\\\
0&\mathbb{i}&0\\\ \end{pmatrix}\begin{pmatrix}\hat{a}_{+1}\\\ \hat{a}_{0}\\\
\hat{a}_{-1}\end{pmatrix}$ (42)
$\hat{J}_{z}=\begin{pmatrix}\hat{a}_{+1}&\hat{a}_{0}&\hat{a}_{-1}\\\
\end{pmatrix}\begin{pmatrix}1&0&0\\\ 0&0&0\\\
0&0&-1\end{pmatrix}\begin{pmatrix}\hat{a}_{+1}\\\ \hat{a}_{0}\\\
\hat{a}_{-1}\end{pmatrix}\ ,$ (43)

where $\mathbb{i}$ is the imaginary unit, $\mathbb{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:

$\displaystyle\hat{\mathbf{J}}^{2}:$
$\displaystyle=\hat{J}_{x}^{2}+\hat{J}_{y}^{2}+\hat{J}_{z}^{2}=$
$\displaystyle=\hat{J}_{z}^{2}+\hat{N}+\hat{N}_{0}+2\hat{N}_{0}(N-\hat{N}_{0})+2((\hat{a}^{\dagger}_{0})^{2}\hat{a}_{+1}\hat{a}_{-1}+\mathrm{h.c.})$
(44)

Alternatively, by identifying the pair annihilation operator
$\hat{\Theta}_{2}=\hat{a}_{0}^{2}-2\hat{a}_{+1}\hat{a}_{-1}$, we can write
$\hat{\textbf{J}}^{2}=N^{2}-\hat{\Theta}^{\dagger}\hat{\Theta}$. As first
observation, we notice that unlike the spin-1/2 case, $\hat{\mathbf{J}}^{2}$
is not diagonal in the occupation basis. Direct computation shows that the
Hamiltonian (39) is equivalent to:

$\hat{H}_{0}(c_{1})=c_{1}\hat{\mathbf{J}}^{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:

$\hat{H}_{0}(c_{2},p)=c_{2}\hat{\bf
J}^{2}+p[\hat{N}_{0}(\hat{N}_{0}-1)+4(\hat{N}_{1}\hat{N}_{-1}+\hat{N}_{-2}\hat{N}_{2})+2(\hat{T}_{0,-2}^{+}\hat{T}_{0,2}^{+}-\hat{T}_{0,-1}^{+}\hat{T}_{0,1}^{+}-2\hat{T}_{1,2}^{+}\hat{T}_{-1,-2}^{+}+(+\leftrightarrow-))]\
,$ (46)

where $c_{2}=(g_{4}-g_{2})/7,p=(7g_{0}-10g_{2}+3g_{4})/35$ and we introduced
the collective operators
$\hat{T}_{s,s^{\prime}}^{+}=\hat{a}^{\dagger}_{s}\hat{a}_{s^{\prime}}+\mathrm{h.c.},\hat{T}_{s,s^{\prime}}^{-}=\mathbb{i}(\hat{a}^{\dagger}_{s}\hat{a}_{s^{\prime}}-\mathrm{h.c.})$.
Using the pair anhiliation operator for spin-2,
$\hat{\Theta}_{2}=\hat{a}_{0}^{2}-2\hat{a}_{1}\hat{a}_{-1}+2\hat{a}_{+2}\hat{a}_{-2}$,
the term coupling $p$ can be interpreted as a pair amplitude
$\hat{\Theta}_{2}^{\dagger}\hat{\Theta}_{2}$. This term commutes with
$\hat{\mathbf{J}}^{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 $\hat{H}=\hat{H}_{0}/N+q\hat{Q}_{z}$, where $q\propto B^{2}$ and
$\hat{Q}_{z}=\sum_{i=1}^{N}\hat{J}_{z}^{2}$. We introduce a scaling to make
the Hamiltonian extensive in $N$. We consider an experiment where all atoms
are initially prepared in the state $\otimes_{i=1}^{N}|0_{z}\rangle$, so that
we may remove the linear Zeeman term, proportional to $\hat{J}_{z}$ (as
$\hat{J}_{z}$ is a symmetry of $\hat{H}$, and acts trivially on the state).

##  Appendix 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:

$\langle\hat{N}_{a,s}\rangle\
,\langle\hat{N}_{a,s}\hat{N}_{a,s^{\prime}}\rangle\ ,\mathrm{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 $\hat{O}_{k}=\sum_{i=1}^{N}\hat{o}_{k}(i)$ with the spin projection
$\hat{J}_{z}$ a symmetry. Since operators acting in different atoms commute,
we notice that the collective operators $\\{\hat{O}_{k}\\}$ obey the same
algebra as their their local counterparts, $\\{\hat{o}_{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\times d$ Hermitian matrices, where
$d=2j+1$. Let $\\{\hat{\lambda}_{\alpha}\\}_{\alpha=1}^{d^{2}-1}$ be a basis
of $d\times d$ Hermitian traceless matrices fulfilling the orthogonality
condition
$\mathrm{Tr}[\hat{\lambda}_{\alpha}\hat{\lambda}_{\beta}]=\delta_{\alpha\beta}$,
for instance the so-called generalized Gell-Mann matrices. The commutation
relations:

$\hat{\lambda}_{\alpha}\hat{\lambda}_{\beta}-\hat{\lambda}_{\beta}\hat{\lambda}_{\alpha}=\mathbb{i}\sum_{\gamma=1}^{d^{2}-1}f_{\alpha\beta\gamma}\hat{\lambda}_{\gamma}\
,$ (48)

are characterized by the structure constants $f_{\alpha\beta\gamma}$. The
structure constants are the same for the generators of the $SU(d)$ group with
the Lie bracket
$[\hat{\lambda}_{\alpha},\hat{\lambda}_{\beta}]:=\hat{\lambda}_{\alpha}\hat{\lambda}_{\beta}-\hat{\lambda}_{\beta}\hat{\lambda}_{\alpha}$.
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
$\\{\hat{\Lambda}_{\alpha}=I[\hat{\lambda}_{\alpha}]\\}_{\alpha=1}^{d^{2}-1}$.
The operator $\hat{\Lambda}_{\alpha}$ is the collective counterpart of
$\hat{\lambda}_{\alpha}$, that is
$\hat{\Lambda}_{\alpha}=\sum_{i=1}^{N}\hat{\lambda}_{\alpha}(i)$.
Consequently, for any other operator $\hat{o}$, the map $I$ distributes
linearly

$I[\hat{o}]=N\mathrm{Tr}[\hat{o}]\frac{\mathbb{I}}{D}+\sum_{\alpha=1}^{d^{2}-1}\mathrm{Tr}[\hat{\lambda}_{\alpha}\hat{o}]\hat{\Lambda}_{\alpha}:=\hat{O}\
,$ (49)

where $D$ is the dimension of the representation, and the identity element
$\mathbb{I}$ is added in order to account for the trace of $\hat{o}$. The
multiplicity $N$ is introduced to consistently transform the identity,
$I[\mathbb{1}]=\sum_{i=1}^{N}\mathbb{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 $\ket{\rm GS}$. Then, the computation of
moments reduces to matrix multiplication matrices and computing expectation
values in the ground state: $\langle\hat{O}_{a}\rangle=\bra{\rm
GS}I[\hat{o}_{a}]\ket{\rm GS}$, $\langle\hat{O}_{a}\hat{O}_{b}\rangle=\bra{\rm
GS}I[\hat{o}_{a}]I[\hat{o}_{b}]\ket{\rm GS}$, etc, for any collective
operators $\hat{O}_{a}$. In this work in particular we probe the GS with
magnetic-sublevel populations along a direction $a$,
$\\{\hat{N}_{a,s}\\}_{s=-j}^{j}$. These are generated from the single-atom
observable $\hat{n}_{a,s}=\prod_{s\neq
s^{\prime}=-j}^{j}(s-\hat{S}_{a})/(s-s^{\prime})$, where $\hat{S}_{a}$ is the
atomic spin projection along direction $a\in\\{x,y,z\\}$, i.e
$\hat{S}_{a}:=\mathbf{a}\cdot\hat{\mathbf{S}}$.

##  Appendix 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 $\hat{X}$ and
$\hat{Y}$, acting on the Hilbert space for $N$ qudits, our goal is find the
possible combinations of “feasible” mean values
$(x=\langle\hat{X}\rangle,y=\langle\hat{Y}\rangle)$, 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 $\langle\hat{Y}\rangle$,
while keeping fixed the mean value $x=\langle\hat{X}\rangle$. 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 $\lambda_{x}$, and finding the ground state
of Hamiltonians of the form:

$\hat{L}(\lambda_{x})=\pm\hat{Y}-\lambda_{x}(\hat{X}-x\mathbb{1})\leavevmode\nobreak\
.$ (50)

Adjusting the multiplier $\lambda_{x}$ ensures that the constraint
$\langle\hat{X}\rangle=x$ is met in the ground state. More generally, the
extremal points $(x,y)$ are found as the expectation values of
$\hat{X},\hat{Y}$ in the ground state of:

$\hat{L}(\lambda_{x},\lambda_{y})=\lambda_{x}\hat{X}+\lambda_{y}\hat{Y}\leavevmode\nobreak\
,$ (51)

varying the parameters $\lambda_{x},\lambda_{y}$ between $-\infty$ and
$+\infty$. 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 $\hat{L}$, the same technique outlined
above for the Bose gas Hamiltonian can be applied since in the cases of
interest here, $\hat{L}$ is a polynomial of collective operators.

##  Appendix 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 $\hat{S}_{a}(i)=\sum_{s=-j}^{j}s\hat{n}_{a,s}(i)$, as well as the
projectors $\hat{n}_{a,0}(s)$ onto the $s=0$ magnetic sublevels, along the
three orthogonal directions $a\in\\{x,y,z\\}$. For completeness, we rewrite
here the EW inequalities (9) and (10) of the main text:

$\mathrm{Var}(\hat{\mathbf{J}})-\frac{1}{N-1}\langle\hat{\mathbf{N}}_{0}(\hat{\mathbf{N}}_{0}-\mathbf{1})\rangle\geq
N[j(j+1)-\beta_{j}]\ ,$ (52)

and:

$(N-1){\rm
Var}(\hat{J}_{z})-N\langle\hat{Q}_{z}\rangle-\langle\hat{J}_{x}^{2}+\hat{J}_{y}^{2}\rangle+N[(N-1)\beta_{j}+j(j+1)]\geq\langle\hat{\mathbf{N}}_{0}(\hat{\mathbf{N}}_{0}-\mathbf{1})\rangle\leavevmode\nobreak\
,$ (53)

where $\mathbf{1}=(1,1,1)$ and $\beta_{j}$ is the corresponding bound in Eq.
(3). We evaluated this bound numerically, and the corresponding values are:

spin $j$ | $1$ | $2$ | $3$ | $4$ | $5$
---|---|---|---|---|---
$\beta$ $\\{\mathbf{S},\mathbf{n}_{0}\\}$ | 1.5 | 4.3 | 9.2 | 16.15 | 25.22
Table 1: Numerical separable bound $\beta$ for the basis of local measurements
$\\{\\{\hat{S}_{a}(i),\hat{n}_{a,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:
$\sum_{a\in\\{x,y,z\\}}(\hat{N}_{a,+1}+\hat{N}_{a,-1})=3N-\sum_{a\in\\{x,y,z\\}}\hat{N}_{a,0}=\sum_{i=1}^{N}[\hat{\mathbf{S}}(i)]^{2}=j(j+1)N=2N$.
Together with the bound $\beta_{j=1}=3/2$, we obtain Eqs. (6), (7) from
respectively Eqs. (52) and (53).

#### Derivation of ${\rm Tr}[C]$.

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

$\displaystyle\mathrm{Tr}[C]=$ $\displaystyle\frac{1}{N(N-1)}\sum_{i\neq
j}\sum_{a\in\\{x,y,z\\}}\left\langle\hat{S}_{a}(i)\hat{S}_{a}(j)+\hat{n}_{a,0}(i)\hat{n}_{a,0}(j)\right\rangle$
(54) $\displaystyle=$
$\displaystyle\frac{1}{N(N-1)}\sum_{a\in\\{x,y,z\\}}\left\langle[\hat{J}_{a}]^{2}+[\hat{N}_{a,0}]^{2}-\sum_{i=1}^{N}[[\hat{S}_{a}(i)]^{2}+\hat{n}_{a,0}(i)]\right\rangle$
(55) $\displaystyle=$
$\displaystyle\frac{1}{N(N-1)}\left(\langle\hat{\mathbf{J}}^{2}+\hat{\mathbf{N}}_{0}(\hat{\mathbf{N}}_{0}-\mathbf{1})\rangle-j(j+1)N\right)$
(56)

where $\hat{\bf J}=(\hat{J}_{x},\hat{J}_{y},\hat{J}_{z})$ is the collective
spin, $\hat{\bf N}_{0}=(\hat{N}_{x,0},\hat{N}_{y,0},\hat{N}_{z,0})$ and ${\bf
1}=(1,1,1)$. We used that $\sum_{a}[\hat{S}_{a}(i)]^{2}=j(j+1)$, and that
$[\hat{n}_{a,0}(i)]^{2}=\hat{n}_{a,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})$.

#### 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}=\begin{pmatrix}1&0\\\
0&0\end{pmatrix}:=P_{j}$ in the local basis described before, leading to

$\displaystyle\mathrm{Tr}[NP(C-\mathbf{m}\otimes\mathbf{m})]=$
$\displaystyle\frac{1}{N-1}\left(\langle\hat{\mathbf{J}}^{2}\rangle-j(j+1)N\right)-\frac{1}{N}\langle\hat{\mathbf{J}}\rangle^{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
$\mathbf{J^{2}}=Nj(Nj+1)$, $J_{z}=0$), the optimal projector is
$P_{x}=P_{y}=0,P_{z}=P_{j}$, from which we find:

$\displaystyle\mathrm{Tr}[NP(C-\mathbf{m}\otimes\mathbf{m})]=$
$\displaystyle\frac{1}{N-1}\langle(\hat{J}_{z})^{2}-\hat{Q}_{z}\rangle-\frac{1}{N}\langle\hat{J}_{z}\rangle^{2}\
.$ (58)

Inserting this expression, together with Eq. (56), into our central inequality
(4), we obtain the announced result Eq. (53).

##  Appendix G Derivation of spin-2 EW, inequality (3) and generalization to
arbitrary integer spins

In this appendix we show how to derive the entanglement witness presented in
Eq. (3), 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:
$\\{\hat{n}_{a,+2}(i),\hat{n}_{a,+1}(i),\hat{n}_{a,-1}(i),\hat{n}_{a,-2}(i)\\}_{a\in\\{x,y,z\\}}$.
This choice leads to the $\mathrm{Tr}(C)$ term:

$\displaystyle\mathrm{Tr}(C)=$
$\displaystyle\frac{1}{N(N-1)}\sum_{a\in\\{x,y,z\\}}\sum_{s\neq
0}\langle[\hat{N}_{a,s}]^{2}-\hat{N}_{a,s}\rangle$ (59)
$\displaystyle\frac{1}{N(N-1)}\sum_{a\in\\{x,y,z\\}}\left(\sum_{s\neq
0}\langle[\hat{N}_{a,s}]^{2}\rangle-(N-\langle\hat{N}_{a,0}\rangle\right)$
(60) $\displaystyle=$
$\displaystyle\frac{1}{N(N-1)}\sum_{a\in\\{x,y,z\\}}\left(\sum_{s\neq
0}\langle[\hat{N}_{a,s}]^{2}\rangle+\langle\hat{N}_{a,0}\rangle\right)-\frac{3}{N-1}$
(61)

The method reveals the optimal projector
$P_{x}=P_{y}=P_{z}=\frac{1}{2}\begin{pmatrix}1&0&0&-1\\\ 0&1&-1&0\\\
0&-1&1&0\\\ -1&0&0&1\end{pmatrix}$. This yields the terms:

$\mathrm{Tr}[NPC]=\frac{1}{2(N-1)}\left(\sum_{a\in\\{x,y,z\\}}\left\langle[\hat{N}_{a,+2}-\hat{N}_{a,-2}]^{2}+[\hat{N}_{a,+1}-\hat{N}_{a,-1}]^{2}+\hat{N}_{a,0}\right\rangle-3N\right)$
(62)
$\mathrm{Tr}[NP(\mathbf{m}\otimes\mathbf{m})]=\frac{1}{2N}\sum_{a\in\\{x,y,z\\}}\left(\langle\hat{N}_{a,+2}-\hat{N}_{a,-2}\rangle^{2}+\langle\hat{N}_{a,+1}-\hat{N}_{a,-1}\rangle^{2}\right)\
.$ (63)

We introduce the notation $\hat{\sigma}_{a,s}=\hat{N}_{a,s}-\hat{N}_{a,-s}$,
the population imbalance between opposed levels along direction $a$. Together
with the notation ${\bf 1}=(1,1,1)$, we can rewrite these equations as:

$2(N-1){\rm Tr}(C)=\frac{2}{N}\sum_{s\neq 0}\langle[\hat{\bf
N}_{s}]^{2}\rangle+\frac{2}{N}\langle\hat{\bf N}_{0}\cdot{\bf 1}\rangle-6$
(64)
$2N(N-1)\mathrm{Tr}[PC]=\langle[\hat{\boldsymbol{\sigma}}_{1}]^{2}+[\hat{\boldsymbol{\sigma}}_{2}]^{2}\rangle+\langle\hat{\bf
N}_{0}\cdot{\bf 1}\rangle-3N$ (65)
$2N(N-1)\mathrm{Tr}[P(\mathbf{m}\otimes\mathbf{m})]=\frac{N-1}{N}[\langle\hat{\boldsymbol{\sigma}}_{1}\rangle^{2}+\langle\hat{\boldsymbol{\sigma}}_{2}\rangle^{2}].$
(66)

Together with the bound $\beta=3/2$ found numerically, we find:

$\displaystyle 2(N-1)\left[N\mathrm{Tr}[P(C-\mathbf{m}\otimes\mathbf{m})]-{\rm
Tr}(C)+\beta\right]=$
$\displaystyle\langle[\hat{\boldsymbol{\sigma}}_{1}]^{2}+[\hat{\boldsymbol{\sigma}}_{2}]^{2}\rangle-\frac{N-1}{N}[\langle\hat{\boldsymbol{\sigma}}_{1}\rangle^{2}+\langle\hat{\boldsymbol{\sigma}}_{2}\rangle^{2}]$
$\displaystyle+\left(1-\frac{2}{N}\right)\langle\hat{\bf N}_{0}\cdot{\bf
1}\rangle-\frac{2}{N}\sum_{s\neq 0}\langle[\hat{\bf
N}_{s}]^{2}\rangle+3\leavevmode\nobreak\ ,$ (67)

namely, Eq. (3).

#### Generalization to arbitrary even-level systems.

A possible generalization to arbitrary integer spin is obtained by taking the
projector to $\sum_{s=1}^{j}(\hat{n}_{a,s}-\hat{n}_{a,-s})$ for
$a\in\\{x,y,z\\}$. The corresponding bound $\beta$ for local measurements
$\\{\\{\hat{n}_{a,s}\\}_{s=-j\neq 0}^{j}\\}_{a\in\\{x,y,z\\}}$ obtained
numerically is summarized in the following table for several spin-$j$.

spin $j$ | 1 | 2 | 3 | 4 | 5
---|---|---|---|---|---
$\beta$ for $\\{\hat{n}\setminus\hat{\mathbf{n}}_{0}\\}$ | 1.25 | 1.5 | 1.5 | 1.3369 | 1.2926
Table 2: Numerical separable bound $\beta$ for the basis of local measurements
$\\{\\{\hat{n}_{a,s}\\}_{s=-j\neq 0}^{j}\\}_{a\in\\{x,y,z\\}}$ for different
spin-$j$.

Notice that for $j=1$, we recover Eq. (6).

##  Appendix 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
$\\{\\{\hat{o}_{k}(i)\\}_{k=1}^{K}\\}_{i=1}^{N}$. 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
$\\{\\{\hat{n}_{a,+1}(i),\hat{n}_{a,-1}(i)\\}_{a=\\{x,y,z\\}}^{K}\\}_{i=1}^{N}$,
where $\hat{n}_{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 $\mathbb{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
$\\{\\{\hat{s}_{a}(i),\hat{n}_{a,0}(i)\\}_{a=\\{x,y,z\\}}^{K}\\}_{i=1}^{N}$.

###  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:

$\left\\{\mathbf{m}_{a}=\frac{1}{N}(\langle\hat{N}_{a,+1}\rangle,\langle\hat{N}_{a,-1}\rangle),\
C_{a}=\frac{1}{N(N-1)}\begin{pmatrix}\langle\hat{N}_{a,+1}(\hat{N}_{a,+1}-1)\rangle&\langle\hat{N}_{a,+1}\hat{N}_{a,-1}\rangle\\\
\langle\hat{N}_{a,-1}\hat{N}_{a,+1}\rangle&\langle\hat{N}_{a,-1}(\hat{N}_{a,-1}-1)\rangle\end{pmatrix}\right\\}_{a\in\\{x,y,z\\}}\
,$ (68)

where $\hat{N}_{a,s}=\sum_{i=1}^{N}\hat{n}_{a,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:

$\mathbf{m}_{x}=\mathbf{m}_{y}=(0.499,0.499),\mathbf{m}_{z}=(0.002,0.002)\
,C_{x}=C_{y}=\begin{pmatrix}0.248&0.251\\\ 0.251&0.248\end{pmatrix}\
,C_{z}=\begin{pmatrix}0.000&0.000\\\ 0.000&0.000\end{pmatrix}\ .$ (69)

One can prepare such state on an unpolarized spin-1 BEC adiabatically by
ramping down the magnetic field to $\sim\sqrt{q}$ from $q/c\gg 1$. Then, one
needs to perform a collective Stern-Gerlach splitting along three orthogonal
spatial directions, as achieved e.g. in ref. [10]. The statistics of (68) can
readily be inferred by counting the number of atoms on each magnetic sublevel
via fluorescence imaging [10].

###  H.3 Detecting quantum entanglement

Now we obtain the matrices
$\\{C_{a}-\mathbf{m}_{a}\otimes\mathbf{m}_{a}:=\tilde{C}_{a}\\}_{a\in\\{x,y,z\\}}$,
which are central to our method. With the values in Eq. (69), we have:

$\tilde{C}_{x}=\tilde{C}_{y}=\begin{pmatrix}-0.139&0.140\\\
0.140&-0.139\end{pmatrix}\ ,\tilde{C}_{z}=\begin{pmatrix}0.000&0.002\\\
0.002&0.000\end{pmatrix}\ .$ (70)

Next, we diagonalize the $\tilde{C}$-matrices to compute the quantity
$\mathcal{W}=\sum_{a\in\\{x,y,z\\}}(\sum_{\lambda<0}\lambda(\tilde{C}_{a})-\mathrm{Tr}(C_{a}))$,
where $\lambda(\tilde{C}_{a})$ are the eigenvalues of $\tilde{C}_{a}$. In this
example, it yields $\mathcal{W}=-1.550$. Comparing this value with the bound
$\beta$ for the above set of local measurements, $\beta=1.25$ (see Table 2 in
Appendix G), we have $\mathcal{W}+\beta=-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\in\\{x,y,,z\\}}$.

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

$P_{x}=P_{y}=P_{z}=\frac{1}{2}\begin{pmatrix}1&-1\\\ -1&1\end{pmatrix}\ .$
(71)

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 $\beta$ 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
$\langle\hat{\mathbf{J}}^{2}\rangle,\langle\hat{\mathbf{J}}\rangle,\langle\hat{\mathbf{N}}_{0}^{2}\rangle$.