# Detecting Bell correlations in multipartite non-Gaussian spin states

Jiajie Guo  State Key Laboratory for Mesoscopic Physics, School of Physics,
Frontiers Science Center for Nano-optoelectronics, $\&$ Collaborative
Innovation Center of Quantum Matter, Peking University, Beijing 100871, China
Jordi Tura  Instituut-Lorentz, Universiteit Leiden, P.O. Box 9506, 2300 RA
Leiden, The Netherlands  Qiongyi He<EMAIL_ADDRESS>State Key Laboratory
for Mesoscopic Physics, School of Physics, Frontiers Science Center for Nano-
optoelectronics, $\&$ Collaborative Innovation Center of Quantum Matter,
Peking University, Beijing 100871, China  Collaborative Innovation Center of
Extreme Optics, Shanxi University, Taiyuan, Shanxi 030006, China  Peking
University Yangtze Delta Institute of Optoelectronics, Nantong 226010,
Jiangsu, China  Matteo Fadel<EMAIL_ADDRESS>Department of Physics, ETH
Zürich, 8093 Zürich, Switzerland

###### Abstract

We expand the toolbox for studying Bell correlations in multipartite systems
by introducing permutationally invariant Bell inequalities (PIBIs) involving
few-body correlators. First, we present around twenty families of PIBIs with
up to three- or four-body correlators, that are valid for arbitrary number of
particles. Compared to known inequalities, these show higher noise robustenss,
or the capability to detect Bell correlations in highly non-Gaussian spin
states. We then focus on finding PIBIs that are of practical experimental
implementation, in the sense that the associated operators require collective
spin measurements along only a few directions. To this end, we formulate this
search problem as a semidefinite program that embeds the constraints required
to look for PIBIs of the desired form.

Some correlations arising from quantum physics cannot be explained within the
paradigm of local realism, and are thus called nonlocal [1]. These are
detected via the violation of a so-called Bell inequality [2], tested in
practice through a Bell experiment [3]. Besides their fundamental interest,
nonlocal correlations are the resource enabling device-independent (DI)
quantum information processing tasks, such as quantum key distribution [4],
randomness amplification [5] or self-testing [6]. Although much research has
focused on few-partite scenarios, mostly bipartite [7, 8], nonlocal
correlations also appear naturally in the multipartite regime [9, 10] and, in
particular, in physically relevant many-body systems [11, 12, 13]. With mild
additional assumptions, multipartite nonlocality can be revealed in
experimentally-practical ways, and take the name of Bell correlations [14,
15].

Detection of Bell correlations is of great interest, as they are related to
quantum critical points [16], metrology [17], open quantum systems [18], and
bosonic systems at finite temperature [19, 20], and provide an avenue to
quantify DI entanglement and Bell correlation depth [21, 22, 23]. However, the
available inequalities are scarce, because a complete characterization is an
intractable task [24]. An approach that finds a good compromise between
expressivity and complexity is to focus on Bell inequalities with particular
symmetries and low-order correlators [11, 12, 25, 26]. In turn, this reduces
the experimental requirements to reveal Bell correlations from them. A
paradigmatic example is the use of two-body, permutationally invariant Bell
inequalities (PIBIs) to detect a class of Gaussian states known as spin-
squeezed states [14, 15].

Despite all this progress, so far only PIBIs with up to two-body correlators
are known, which poses a fundamental limit on their applicability. It is thus
of great interest to find PIBIs involving higher-order moments of physical
observables, such that Bell correlations can be detected in larger classes of
states and with higher noise tolerance. In particular, these tools would
enable the study of Bell correlations in non-Gaussian states [27], which
cannot be characterized by only second moments.

In this work we present around twenty new PIBIs involving three- and four-body
correlators, and illustrate that compared to known PIBIs they provide an
advantage in terms of noise robustness and sensitivity to non-Gaussian states.
Moreover, we provide a general framework to derive new PIBIs with high-order
correlators under the constraint of being experimentally practical, in the
sense that they can be tested by performing collective spin measurements along
only a few directions. This is based on a semidefinite program (SDP) that
allows us to find Bell correlation witnesses of a desired ansatz.

Preliminaries.– We consider the multipartite Bell experiment in which N
observers, labeled as $i=1,2,\cdots,N$, perform one of the two measurements
$M_{0}^{(i)},M_{1}^{(i)}$, on their part of the system, and obtain one of the
two possible outcomes $\pm 1$. Correlations among parties are characterized by
the correlators $\langle M_{j_{1}}^{(i_{1})}\cdots
M_{j_{k}}^{(i_{k})}\rangle$. However, to reduce the complexity of
characterising multipartite correlations, we restrict ourselves to
permutationally invariant (PI) observables, namely

$\displaystyle\mathcal{S}_{j_{1}\cdots j_{k}}$
$\displaystyle=\sum_{\stackrel{{\scriptstyle 1\leq i_{1}<\cdots<i_{k}\leq
N}}{{\text{+ perm.}}}}\langle M_{j_{1}}^{(i_{1})}\cdots
M_{j_{k}}^{(i_{k})}\rangle\;.$ (1)

Here, $k=1,...,K$ indicates the order of the PI correlator, and $j_{l}=0,1$ is
the measurement setting for $l=1,\ldots,k$.

By considering $N$ parties and at most $K$th-order PI correlators, the set of
classical correlations form a polytope $\mathbb{P}_{N,K}^{S}$.The vertices of
this polytope are identified with the correlations originating from
deterministic local hidden variable models (LHVM), for which it holds $\langle
M_{j_{1}}^{(i_{1})}\cdots M_{j_{k}}^{(i_{k})}\rangle=\langle
M_{j_{1}}^{(i_{1})}\rangle\cdots\langle M_{j_{k}}^{(i_{k})}\rangle$ and
$\langle M_{j}^{(i)}\rangle=\pm 1$ for all $i,j$. From these vertices, it is
possible to derive a dual description of $\mathbb{P}_{N,K}^{S}$ in terms of
the linear inequalities (i.e. PIBIs) defining its facets. Measurement
statistics lying outside $\mathbb{P}_{N,K}^{S}$ indicate the presence of Bell
correlations. For details about this framework we refer the reader to Refs.
[11, 12].

In a scenario with $N<20$ and $K\lesssim 3$ it is relatively simple to list
all vertices of $\mathbb{P}_{N,K}^{S}$ and to obtain from them the full list
of facet inequalities. For larger $N$ and $K$, however, this approach is
unfeasible. Therefore, to find PIBIs that allow the detection of Bell
correlations in many-body systems, we might rely on the following method. For
a given $K$, we characterize $\mathbb{P}_{N,K}^{S}$ for a few small values of
$N$. Then, we look for patterns in the inequalities that appear, and use them
to conjecture families of PIBIs valid for arbitrary $N$. Finally, we prove the
conjectured families by demonstrating that they cannot be violated by LHVM.

Until now, only a couple of PIBIs valid for arbitrary $N$ are known, and only
for $K=2$ [11, 28]. Of particular relevance is the inequality

$I_{2}\equiv-2\mathcal{S}_{0}+\frac{1}{2}\mathcal{S}_{00}-\mathcal{S}_{01}+\frac{1}{2}\mathcal{S}_{11}+2N\geq
0\;,$ (2)

which enabled the experimental detection of Bell correlations in spin-squeezed
BECs [14] and cold atomic ensembles [15].

Third-order Bell inequalities.– We start with considering the case $K=3$.
Remarkably, by computing the polytope $\mathbb{P}_{N,3}^{S}$ for all $N<12$,
we were able to identify around twenty families of Bell inequalities (See
Supplementary Sec. V). As we will see, an interesting family among these is

$\displaystyle
I_{3}\equiv-12(N-1)\mathcal{S}_{0}-12(N-1)\mathcal{S}_{1}+3(N-2)\mathcal{S}_{00}+6N\mathcal{S}_{01}+$
$\displaystyle\;3(N-2)\mathcal{S}_{11}-2\mathcal{S}_{000}-3\mathcal{S}_{001}+\mathcal{S}_{111}+12N(N-1)\geq
0$ (3)

which we have proven to be valid for all $N$ in Supplementary Sec. I.

As all the inequalities we consider here involve two measurement settings and
two outcomes, Jordan’s lemma guarantees that their maximum quantum violation
can be achieved by qubit measurements [29, 30]. To this end, we identify
$M_{j}^{(i)}=\vec{u}_{j}\cdot\vec{\sigma}^{(i)}\;(j=0,1)$ for the $i$th party,
where $\vec{\sigma}$ is the vector of Pauli matrices. Even if not necessarily
optimal nor required by the inequality, we assume that the same pair of
observables is chosen by all parties, i.e. $M_{j}^{(i)}=M_{j}$. Since local
rotations are irrelevant here, we can simplify further our discussion by
choosing $M_{0}=\sigma_{z}$ and
$M_{1}(\theta)=\sin(\theta)\,\sigma_{x}+\cos(\theta)\,\sigma_{z}$, where
$\theta\in[0,\pi]$. With this definition, the correlators Eq. (1) can be
written as operators

$\displaystyle\hat{\mathcal{S}}_{j_{1}\cdots j_{k}}(\theta)$
$\displaystyle=\sum_{\stackrel{{\scriptstyle 1\leq i_{1}<\cdots<i_{k}\leq
N}}{{\text{+ perm.}}}}M_{j_{1}}^{(i_{1})}(\theta)\otimes\cdots\otimes
M_{j_{k}}^{(i_{k})}(\theta)\;,$ (4)

so that the Bell inequalities (2) and (Detecting Bell correlations in
multipartite non-Gaussian spin states) can now be understood as Bell operators
$\hat{I}_{2}(\theta)$ and $\hat{I}_{3}(\theta)$, respectively.

Given a Bell operator, we search for the optimal measurements and states that
maximizes the quantum violation relative to the classical bound, i.e.
$Q_{V}^{N}/\beta_{C}^{N}$, where the classical bound $\beta_{C}^{N}$ is the
constant term appearing in the inequalities. Bell nonlocality is detected for
a state $|\psi\rangle$ if the Bell operator yields a negative expectation
value $\langle\psi|\hat{I}_{K}|\psi\rangle<0$. As the classical bound is a
constant, the maximum quantum violation $Q_{V}^{N}$ can thus be identified
with the minimum eigenvalue of Bell operator

$\displaystyle
Q_{V}^{N}=\min_{\theta}\min_{|\psi\rangle}\langle\psi|\hat{I}_{K}|\psi\rangle=\min_{\theta}\lambda_{min}(\hat{I}_{K})\;,$
(5)

and the associated eigenvector is the state that maximally violates $I_{K}$.

In the $N$-qubit Hillbert space, the dimension of Bell operator scales
exponentially with $N$, making challenging to solve the eigenvalue problem Eq.
(5) for large $N$. Fortunately, since the correlators $\mathcal{S}$ have
permutation symmetry, it is possible to introduce a symmetry-adapted basis in
which to express the Bell operator, such that it block-diagonalizes due to
Schur-Weyl duality [31, 32]. We can then focus the search of nonlocality onto
the fully symmetric block $\hat{I}_{K}^{\text{Sym}}$ of size
$(N+1)\times(N+1)$, and the maximum quantum violation of $\hat{I}_{K}$ can
then be obtained from the lowest eigenvalue of $\hat{I}_{K}^{\text{Sym}}$
[12]. In Fig.1 we show $Q_{V}^{N}/\beta_{C}$ for $I_{2}$ and $I_{3}$, as a
function of $N$. It is evident the significantly better scaling for the
higher-order Bell inequalities $I_{3}$ compared to $I_{2}$. In the limit
$N\rightarrow\infty$, it is possible to show through a variational calculation
[12] that the relative violation of $I_{3}$ tends to
$-2\sqrt{3}/9\approx-0.3849$, which is larger than the value $-1/4$ obtained
for $I_{2}$. This indicates a higher noise robustness, as well as the
possibility to detect Bell correlations in a larger class of states, as we
will see with an example in the next paragraph. A similar analysis for the
quantum violation of some other third-order Bell inequalities $I_{3}$ and a
fourth-order Bell inequality $I_{4}$ are given in Supplementary Sec. V, VI.

Figure 1: Maximum relative quantum violation $Q_{V}^{N}/\beta_{C}^{N}$ for the
3rd-order Bell inequality $I_{3}$ and 2nd-order Bell inequality $I_{2}$, as a
function of the number of parties $N$. The two horizontal dashed lines
indicate the asymptotic violation for $N\rightarrow\infty$, which for $I_{2}$
is $-1/4$, and for $I_{3}$ is $-2\sqrt{3}/9\approx-0.3849$.

Bell correlations in spin-squeezed states.– We now show that $I_{3}$ allows us
to detect Bell correlations in many-body spin states of experimental
relevance. To illustrate this, let us consider the states that can be prepared
through the paradigmatic one-axis twisting (OAT) Hamiltonian
$H=\hbar\chi\hat{S}_{z}^{2}$ [33]. The evolution of an initial coherent spin
state along the $x$-axis for a time $t$ can be parametrized through the
(adimensional) interaction strength $\mu=2\chi t$, and reads
$|\Phi(\mu)\rangle=2^{-N/2}\sum_{k=0}^{N}\sqrt{\binom{N}{k}}e^{-i(N/2-k)^{2}\mu/2}|k\rangle$,
where $|k\rangle$ represents the $N$-particle Dicke state with $k$
excitations.

To investigate Bell correlations in the OAT states $|\Phi(\mu)\rangle$, we
compute $I_{2}$ and $I_{3}$ as a function of $\mu$. For this, we need to
minimize the associated Bell operators over the measurement directions $M_{j}$
at every $\mu$, which in the case of a given state have to be both
parametrized over the full sphere as $M_{j}=\vec{u}_{j}\cdot\vec{\sigma}$,
with
$\vec{u}_{j}=(\cos(\phi_{j})\sin(\theta_{j}),\sin(\phi_{j})\sin(\theta_{j}),\cos(\theta_{j}))$.
This yields Bell operators
$\hat{I}_{K}(\phi_{0},\theta_{0},\phi_{1},\theta_{1})$ that are now functions
of four angles.

Concretely, we express the Bell operators in the fully-symmetric subspace in
terms of collective spin operators, make use of the analytic results for their
expectation values for spin OAT states (see Supplementary Sec. II, III), and
then minimize the lowest eigenvalue of the operator over the measurement
directions. The violations we obtain are shown in Fig. 2 for $N=50$, where we
can observe that $I_{3}$ outperforms $I_{2}$ by reaching larger relative
violation $Q_{V}^{N}/\beta_{C}^{N}$ as well as detecting Bell correlations
over a wider squeezing range, and thus for a larger class of states. To
investigate the noise robustness, we consider OAT states mixed with white
noise as
$\rho(\eta,\mu)=\eta\ket{\Phi(\mu)}\bra{\Phi(\mu)}+(1-\eta)\mathbf{I}/(N+1)$.
In Fig. 2 we plot the minimum $\eta$ for observing a PIBI violation, and show
that high-order inequalities detect Bell correlations with higher noise
tolerance.

Figure 2: a): Relative quantum violation of the PIBIs $I_{2}$ (blue) and
$I_{3}$ (orange) for $N=50$ spin OAT states $|\psi(\mu)\rangle$ as a function
of $\mu$. An advantage over $I_{2}$ can also be found for $I_{3}^{\text{SDP}}$
(red dashed), which requires to measure only one third moment of the
collective spin. b): For mixed states
$\rho(\eta,\mu)=\eta\ket{\Phi(\mu)}\bra{\Phi(\mu)}+(1-\eta)\mathbf{I}/(N+1)$,
the minimum purity $\eta$ required to violate each PIBI.

Finding practical high-order Bell inequalities.– We have shown, taking $I_{3}$
as an example, that high-order Bell inequalities can outperform their low-
order counterparts and allow us to robustly detect Bell correlations in states
that are routinely investigated experimentally. However, we note that the
associated Bell operator often requires to measure high-order moments of the
collective spin along several directions, because we have e.g. (see
Supplementary Sec. II)

$\mathcal{S}_{001}=\frac{8}{3}\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\rangle+...\;.$
(6)

This can pose practical challenges, as estimating many high-order moments
requires to collect large measurement statistics. For this reason, we would
like to find inequalities with coefficients satisfying additional constraints,
such that the associated operator involves e.g. only _one_ third-order moment.
To see the form of such constraint, note that for some measurement direction
$\vec{a}=\alpha\vec{m}+\beta\vec{n}$ we have (see Supplementary Sec. IV)

$\langle\hat{S}^{3}_{\vec{a}}\rangle=\frac{\beta^{3}}{8}\mathcal{S}_{000}+\frac{\alpha^{3}}{8}\mathcal{S}_{111}+\frac{3\alpha\beta^{2}}{8}\mathcal{S}_{001}+\frac{3\alpha^{2}\beta}{8}\mathcal{S}_{011}+f(\mathcal{S}_{0},\mathcal{S}_{1})\;,$
(7)

where $f(\mathcal{S}_{0},\mathcal{S}_{1})$ is a linear function of one-body
correlators only, and $\alpha,\beta\in\mathbb{R}$,
$|\alpha^{2}|+|\beta^{2}|=1$. Therefore, only inequalities whose coefficients
for the third-order correlators are following the pattern of Eq. (7) will
result in a Bell operator involving only one third-order moment of the
collective spin.

Note here that imposing such constraints as a further projection of the local
polytope is not trivial, because of the high nonlinearity in $\alpha$ and
$\beta$ of the coefficients multiplying the correlators. One obtains a
different polytope projection for each pair $(\alpha,\beta)$, which
dramatically increases the optimization complexity, as for each projection new
families of Bell inequalities need to be found. Moreover, such approach does
not guarantee to find tight inequalities, as these might correspond to tilts
of facets.

We thus propose a way to circumvent these difficulties, by developing a method
to find practical Bell operators that is based on a hierarchy of semidefinite
programs (SDPs) that search for inequalities of a particular form. For the
sake of brevity and clarity, in the following we present our method in brief,
and show a concrete example for $K=3$. For a more detailed and general
formulation we refer the reader to Supplementary Sec. VII.

To illustrate our idea, let us start recalling that for deterministic LHVM the
PI correlators Eq. (1) can be written as polynomials in four non-negative
integers $(a,b,c,d)$ such that $a+b+c+d=N$ (see Supplementary Sec. I). This
allows us to parametrize the vertices of $\mathbb{P}_{N,K}^{S}$ with integer
partitions of $N$, but it also implies that considering $(a,b,c,d)$ to be non-
negative reals gives us a outer approximation of $\mathbb{P}_{N,K}^{S}$ in
terms of a semialgebraic set. In the space of PI correlators, the latter is
specified by: (i) the set of polynomial equalities $f_{i}(\vec{\cal S})=0$
expressing constraints among correlators, e.g. ${\cal S}_{000}={\cal
S}_{0}^{3}-(3(\mathcal{S}_{0}^{2}-\mathcal{S}_{00})-2){\cal S}_{0}$, and (ii)
the set of four polynomial inequalities $g_{i}(\vec{\cal S})\geq 0$ expressing
$(a,b,c,d)\geq 0$ through correlators, e.g. $a\geq 0$ implies $g_{1}(\vec{\cal
S})=(\mathcal{S}_{0}^{2}-\mathcal{S}_{00})+\mathcal{S}_{0}+\mathcal{S}_{1}+(\mathcal{S}_{0}\mathcal{S}_{1}-\mathcal{S}_{01})\geq
0$.

Having at hand an outer approximation of $\mathbb{P}_{N,K}^{S}$ in terms of a
semialgebraic set, we can apply known techniques that use SDP hierarchies to
test membership to the convex hull of such set [34, 35, 36]. The first step
consists in defining a basis vector, e.g.
$\vec{b}=(1,\mathcal{S}_{0},\cdots,\mathcal{S}_{111})^{T}$, from which the
moment matrix
$\tilde{\Gamma}=\displaystyle\oplus_{i=1}^{4}g_{i}\vec{b}\cdot\vec{b}^{T}$ is
constructed. Then, to check whether the point $\vec{\mathcal{S}}^{*}$ is
outside the convex hull approximating $\mathbb{P}_{N,K}^{S}$ we can linearize
$\tilde{\Gamma}$ and write the SDP

$\begin{array}[]{cccc}\max_{\tilde{\Gamma}}&\lambda&\\\
\mathrm{s.t.}&\tilde{\Gamma}&\succeq&0\\\ &\tilde{\Gamma}_{00}&=&1\\\
&\tilde{\Gamma}_{0i}&=&\lambda(\vec{\mathcal{S}}^{*})_{i}\\\
&\tilde{\Gamma}_{ij}&=&p(\tilde{\Gamma})\;,\end{array}$ (8)

where $p(\tilde{\Gamma})$ is the function expressing the constraints between
the entries of $\tilde{\Gamma}$. If SDP (8) returns $\lambda<1$ we must
conclude that $\vec{\mathcal{S}}^{*}_{3}$ lies outside the convex hull outer
approximating $\mathbb{P}_{N,K}^{S}$, and therefore that this point cannot be
described by a LHVM. In this case, the SDP dual to (8) gives us a certificate
of nonlocality by providing a PIBI that is violated by
$\vec{\mathcal{S}}^{*}_{3}$.

To ensure that the dual to (8) provides us with PIBIs that are experimentally
practical, we now modify SDP (8) by adding the additional constraints required
to obtain Bell operators of the desired form. Note however that this is
nontrivial, as such constraints are highly nonlinear (cf. Eq. (7)). For
example, asking for a Bell operator involving only one third moment of the
collective spin requires to impose the constraints

$\displaystyle(\tilde{\Gamma}_{01},\cdots,\tilde{\Gamma}_{05})$
$\displaystyle=\lambda(\mathcal{S}_{0}^{*},\mathcal{S}_{1}^{*},\cdots,\mathcal{S}_{11}^{*})$
(9) $\displaystyle y$
$\displaystyle=\lambda(\beta^{3}\mathcal{S}_{000}^{*}+3\alpha\beta^{2}\mathcal{S}_{001}^{*}+3\alpha^{2}\beta\mathcal{S}_{011}^{*}+\alpha^{3}\mathcal{S}_{111}^{*})\;,$

where
$y=(\beta^{3}\tilde{\Gamma}_{06}+3\alpha\beta^{2}\tilde{\Gamma}_{07}+3\alpha^{2}\beta\tilde{\Gamma}_{08}+\alpha^{3}\tilde{\Gamma}_{09})$.
Therefore, we run SDP (8) with Eq. (9) as an optimization problem over
$(\alpha,\beta)$, in order to find $\min_{\alpha,\beta}(\lambda)$.

Example with OAT states.– Let us go back to the problem of detecting Bell
correlations in OAT spin states. We aim to find a third-order PIBI with
coefficients such that the resulting Bell operator requires the measurement of
only one third moment of the collective spin.

First, we specify the target state $|\Phi(\mu)\rangle$ where to detect Bell
correlations, e.g. by choosing $N=50$ and $\mu=0.2$. Then, we find a pair of
measurement axes $\vec{u}_{1,2}$ such that the resulting list of correlators
$\vec{\mathcal{S}}^{*}_{3}=(\mathcal{S}_{0}^{*},\mathcal{S}_{1}^{*},\mathcal{S}_{00}^{*},\cdots,\mathcal{S}_{111}^{*})$
shows Bell correlations. This step can be implemented by using SDP (8) to run
the optimization problem $\min_{\vec{u}_{1},\vec{u}_{2}}(\lambda)$. At this
point, note that the dual to SDP (8) will provide us a PIBI that is violated
by $\vec{\mathcal{S}}^{*}_{3}$, but that is in general of difficult
experimental implementation. For this reason, we modify SDP (8) to include the
constraints Eq. (9) and run the optimization problem
$\min_{\alpha,\beta}(\lambda)$. The dual SDP gives now a PIBI in the form

$\displaystyle
I_{3}^{\text{SDP}}=c_{0}+c_{1}\mathcal{S}_{0}+c_{2}\mathcal{S}_{1}+\cdots+c_{6}\left(\beta^{3}\mathcal{S}_{000}+...\right)\geq
0\;,$ (10)

where $\vec{c}=(c_{0},\cdots,c_{5},c_{6})$ are the variables dual to
$(\tilde{\Gamma}_{00},\cdots,\tilde{\Gamma}_{05},y)$ respectively, and whose
associated Bell operator involves only one third-order moment of the
collective spin. For our target state, we obtain $\alpha/\beta=41/59$ and
$\vec{c}=(1,-0.0055,-0.0141,0.0046,0.0099,0.0051,-56.1412)$. The resulting
$I_{3}^{\text{SDP}}$ has a worse noise tolerance than $I_{3}$, but it still
outperforms $I_{2}$, see Fig. 2. This could be improved by searching for a
different inequality for each $\mu$, or by allowing the measurement of two
third-order moments.

Fourth-order Bell inequalities.– Following the ideas presented so far, we can
now guide our search for experimentally practical high-order PIBIs. For $K=4$
we find the inequality

$\displaystyle I_{4}\equiv
24(N-1)\mathcal{S}_{00}+48(N-1)\mathcal{S}_{01}+24(N-3)\mathcal{S}_{11}+$ (11)
$\displaystyle+\mathcal{S}_{0000}+4\mathcal{S}_{0001}+6\mathcal{S}_{0011}+4\mathcal{S}_{0111}+\mathcal{S}_{1111}+48N(N-1)\geq
0$

which involves two- and four-body correlators. The associate Bell operator
reads

$\displaystyle\hat{I}_{4}$
$\displaystyle=64\hat{S}^{4}_{\vec{a}}-192\hat{S}^{2}_{\vec{m}}+32\,g(\vec{a},\vec{m})\,\hat{S}^{2}_{\vec{a}}+h(\vec{a},\vec{m})$
(12)

where $g$ and $h$ are scalar functions of $\vec{a}\cdot\vec{m}$ (see
Supplementary Sec. VI). Remarkably, note that Eq. (12) involves measurements
of the collective spin operator along two directions only. Diagonalising
$\hat{I}_{4}$ according to Eq. (5) we further conclude that: (i) its maximum
relative quantum violation increases with $N$, (ii) the states maximally
violating inequality (11) are highly non-Gaussian, and resembling a
superposition of OAT states, see Fig. 3. The latter states do not violate
$I_{2}$, as it is also expected from the fact that they have zero polarisation
(see Supplementary Sec. VIII).

Figure 3: For $N=50$, Wigner function of the eigenstate corresponding to the
minimum eigenvalue of Bell operator Eq. (12). Red dots indicate the optimal
measurement direction $\hat{S}_{\vec{n}},\hat{S}_{\vec{m}}$, for violating
$I_{4}$.

Conclusions.– We addressed the problem of finding multipartite PIBIs involving
correlators of order higher than second, and that are of practical
experimental implementation in the sense that the associated operators require
collective measurements along only a few directions. We propose about twenty
new PIBIs valid for arbitrary $N$ that involve up to three-body correlators,
and one that involves two- and four-body correlators. From a systematic
analysis, we conclude that in general these inequalities indeed outperform the
currently known PIBIs, since they show higher noise tolerance and the ability
to detect Bell correlations in highly non-Gaussian states. In general, for a
PIBI to be experimentally practical, we note that the coefficients of the
correlators must satisfy some (nonlinear) constraints. We find that these can
be imposed a priori, and formulate a SDP that looks for PIBIs resulting in
Bell operators of the desired form (e.g. involving only one third-moment). Our
results can pave the way to studying Bell correlations in non-Gaussian spin
states, and to use generalized spin-squeezing parameters as Bell correlation
witnesses.

Acknowledgments.– This work is supported by the National Natural Science
Foundation of China (Grants No. 12125402 and No. 11975026). JG acknowledges
financial support from the China Scholarship Council (Grant No. 202106010192).
QH acknowledges the Beijing Natural Science Foundation (Z190005) and the Key
R$\&$D Program of Guangdong Province (Grant No. 2018B030329001). MF was
supported by The Branco Weiss Fellowship – Society in Science, administered by
the ETH Zürich.

## References

  * Brunner  _et al._ [2014] N. Brunner, D. Cavalcanti, S. Pironio, V. Scarani, and S. Wehner, Rev. Mod. Phys. 86, 419 (2014).
  * Bell [1964] J. S. Bell, Physics 1, 195 (1964).
  * Hensen  _et al._ [2015] B. Hensen, H. Bernien, A. E. Dréau, A. Reiserer, N. Kalb, M. S. Blok, J. Ruitenberg, R. F. L. Vermeulen, R. N. Schouten, C. Abellán, W. Amaya, V. Pruneri, M. W. Mitchell, M. Markham, D. J. Twitchen, D. Elkouss, S. Wehner, T. H. Taminiau, and R. Hanson, Nature 526, 682 (2015).
  * Acín  _et al._ [2007] A. Acín, N. Brunner, N. Gisin, S. Massar, S. Pironio, and V. Scarani, Phys. Rev. Lett. 98, 230501 (2007).
  * Colbeck and Renner [2012] R. Colbeck and R. Renner, Nature Physics 8, 450 (2012).
  * Šupić and Bowles [2020] I. Šupić and J. Bowles, Quantum 4, 337 (2020).
  * Pitowsky and Svozil [2001] I. Pitowsky and K. Svozil, Phys. Rev. A 64, 014102 (2001).
  * Rosset  _et al._ [2014] D. Rosset, J.-D. Bancal, and N. Gisin, J. Phys. A: Math. Theor. 47, 424022 (2014).
  * Mermin [1990] N. D. Mermin, Phys. Rev. Lett. 65, 1838 (1990).
  * Tóth  _et al._ [2006] G. Tóth, O. Gühne, and H. J. Briegel, Phys. Rev. A 73, 022303 (2006).
  * Tura  _et al._ [2014a] J. Tura, R. Augusiak, A. B. Sainz, T. Vértesi, M. Lewenstein, and A. Acín, Science 344, 1256 (2014a).
  * Tura  _et al._ [2015] J. Tura, R. Augusiak, A. Sainz, B. Lücke, C. Klempt, M. Lewenstein, and A. Acín, Annals of Physics 362, 370 (2015).
  * Kitzinger  _et al._ [2021] J. Kitzinger, X. Meng, M. Fadel, V. Ivannikov, K. Nemoto, W. J. Munro, and T. Byrnes, Phys. Rev. A 104, 043323 (2021).
  * Schmied  _et al._ [2016] R. Schmied, J.-D. Bancal, B. Allard, M. Fadel, V. Scarani, P. Treutlein, and N. Sangouard, Science 352, 441 (2016).
  * Engelsen  _et al._ [2017] N. J. Engelsen, R. Krishnakumar, O. Hosten, and M. A. Kasevich, Phys. Rev. Lett. 118, 140401 (2017).
  * Piga  _et al._ [2019] A. Piga, A. Aloy, M. Lewenstein, and I. Frérot, Phys. Rev. Lett. 123, 170604 (2019).
  * Fröwis  _et al._ [2019] F. Fröwis, M. Fadel, P. Treutlein, N. Gisin, and N. Brunner, Phys. Rev. A 99, 040101 (2019).
  * Marconi  _et al._ [2022] C. Marconi, A. Riera-Campeny, A. Sanpera, and A. Aloy, Phys. Rev. A 105, L060201 (2022).
  * Fadel and Tura [2018] M. Fadel and J. Tura, Quantum 2, 107 (2018).
  * Niezgoda  _et al._ [2019] A. Niezgoda, J. Chwedeńczuk, L. Pezzé, and A. Smerzi, Phys. Rev. A 99, 062115 (2019).
  * Aloy  _et al._ [2019] A. Aloy, J. Tura, F. Baccari, A. Acín, M. Lewenstein, and R. Augusiak, Phys. Rev. Lett. 123, 100507 (2019).
  * Tura  _et al._ [2019] J. Tura, A. Aloy, F. Baccari, A. Acín, M. Lewenstein, and R. Augusiak, Phys. Rev. A 100, 032307 (2019).
  * Baccari  _et al._ [2019] F. Baccari, J. Tura, M. Fadel, A. Aloy, J.-D. Bancal, N. Sangouard, M. Lewenstein, A. Acín, and R. Augusiak, Phys. Rev. A 100, 022121 (2019).
  * Babai  _et al._ [1991] L. Babai, L. Fortnow, and C. Lund, Computational Complexity 1, 3 (1991).
  * Tura  _et al._ [2014b] J. Tura, A. B. Sainz, T. Vértesi, A. Acín, M. Lewenstein, and R. Augusiak, J. Phys. A: Math. Theor. 47, 424024 (2014b).
  * Tura  _et al._ [2017] J. Tura, G. De las Cuevas, R. Augusiak, M. Lewenstein, A. Acín, and J. I. Cirac, Phys. Rev. X 7, 021005 (2017).
  * Walschaers [2021] M. Walschaers, PRX Quantum 2, 030204 (2021).
  * Wagner  _et al._ [2017] S. Wagner, R. Schmied, M. Fadel, P. Treutlein, N. Sangouard, and J.-D. Bancal, Phys. Rev. Lett. 119, 170403 (2017).
  * [29] B. Toner and F. Verstraete, arXiv:quant-ph/0611001 .
  * Jordan [1870] C. Jordan,  _Traité des substitutions et des équations algébriques_ (Gauthier-Villars, 1870).
  * William Fulton [1999] J. H. William Fulton,  _Representation Theory_ (Springer New York, 1999).
  * Moroder  _et al._ [2012] T. Moroder, P. Hyllus, G. Tóth, C. Schwemmer, A. Niggebaum, S. Gaile, O. Gühne, and H. Weinfurter, New J. Phys. 14, 105001 (2012).
  * Kitagawa and Ueda [1993] M. Kitagawa and M. Ueda, Phys. Rev. A 47, 5138 (1993).
  * Grigoriy Blekherman [2013] R. T. Grigoriy Blekherman, Pablo A. Parrilo, ed.,  _Semidefinite Optimization and Convex Algebraic Geometry_ (SOC FOR INDUSTRIAL & APPLIED M, 2013).
  * Gouveia  _et al._ [2010] J. Gouveia, P. A. Parrilo, and R. R. Thomas, SIAM J. Optim. 20, 2097 (2010).
  * Gouveia and Thomas [2012] J. Gouveia and R. R. Thomas, in  _Semidefinite Optimization and Convex Algebraic Geometry_ (Society for Industrial and Applied Mathematics, 2012) pp. 293–340.
  * Fine [1982] A. Fine, Phys. Rev. Lett. 48, 291 (1982).
  * Fadel and Tura [2017] M. Fadel and J. Tura, Phys. Rev. Lett. 119, 230402 (2017).
  * Lasserre [2001] J. B. Lasserre, SIAM J. Optim. 11, 796 (2001).
  * Walschaers  _et al._ [2018] M. Walschaers, S. Sarkar, V. Parigi, and N. Treps, Phys. Rev. Lett. 121, 220501 (2018).
  * Ra  _et al._ [2020] Y.-S. Ra, A. Dufour, M. Walschaers, C. Jacquard, T. Michel, C. Fabre, and N. Treps, Nat. Phys. 16, 144 (2020).
  * Davis  _et al._ [2021] J. Davis, M. Kumari, R. B. Mann, and S. Ghose, Phys. Rev. Res. 3, 033134 (2021).
  * Kenfack and Życzkowski [2004] A. Kenfack and K. Życzkowski, J. Opt. B: Quantum Semiclass. Opt. 6, 396 (2004).
  * Mari and Eisert [2012] A. Mari and J. Eisert, Phys. Rev. Lett. 109, 230503 (2012).

##  I Supplemental material for “Detecting Bell correlations in multipartite
non-Gaussian spin states”

##  II I. Expressions of $\mathcal{S}$ for local deterministic strategies

For our scenario, in which parties have two measurement settings
$\\{M_{0},M_{1}\\}$ and two possible outcomes $\\{\pm 1\\}$, only four local
deterministic strategies (LDS) are possible, corresponding to the combinations
$\langle M_{0(1)}\rangle=\pm 1$. Since PIBIs are invariant under permutation
of the parties, it is not important which party follows which strategy, but
only the number of parties following each of the strategies. For this reason,
we define the four non-negative integer variables

$\displaystyle a=\\#\\{i\in\\{1,\cdots,N\\}\,|\,\langle
M_{0}^{(i)}\rangle=1,\langle M_{1}^{(i)}\rangle=1\\},$ $\displaystyle
b=\\#\\{i\in\\{1,\cdots,N\\}\,|\,\langle M_{0}^{(i)}\rangle=1,\langle
M_{1}^{(i)}\rangle=-1\\},$ $\displaystyle
c=\\#\\{i\in\\{1,\cdots,N\\}\,|\,\langle M_{0}^{(i)}\rangle=-1,\langle
M_{1}^{(i)}\rangle=1\\},$ $\displaystyle
d=\\#\\{i\in\\{1,\cdots,N\\}\,|\,\langle M_{0}^{(i)}\rangle=-1,\langle
M_{1}^{(i)}\rangle=-1\\},$ (13)

counting the number of parties using each one of the four strategies, such
that $a+b+c+d=N$.

Using the fact that for LDS correlators factorise as

$\displaystyle\langle M^{(i)}_{j_{1}}\cdots
M^{(j)}_{j_{k}}\rangle_{\text{LDS}}=\langle
M^{(i)}_{j_{1}}\rangle\cdot\dots\cdot\langle M^{(j)}_{j_{k}}\rangle\;,$ (14)

we can express PI correlators $\mathcal{S}$ as polynomials in $\\{a,b,c,d\\}$.
For example, one-body correlators read

$\displaystyle\mathcal{S}_{k}=\sum_{i=1}^{N}\langle M^{(i)}_{k}\rangle\;,$
(15)

such that we have the polynomials

$\displaystyle\mathcal{S}_{0}$ $\displaystyle=a+b-c-d\;,$ (16)
$\displaystyle\mathcal{S}_{1}$ $\displaystyle=a-b+c-d\;.$ (17)

For the two-body correlators we can write

$\displaystyle\mathcal{S}_{kl}$ $\displaystyle=\sum_{i,j=1(i\neq
j)}^{N}\langle M^{(i)}_{k}M^{(j)}_{l}\rangle=\sum_{i,j=1}^{N}\langle
M^{(i)}_{k}M^{(j)}_{l}\rangle-\sum_{i=1}^{N}\langle
M^{(i)}_{k}M^{(i)}_{l}\rangle$ $\displaystyle=\sum_{i=1}^{N}\langle
M^{(i)}_{k}\rangle\sum_{j=1}^{N}\langle
M^{(j)}_{l}\rangle-\sum_{i=1}^{N}\langle M^{(i)}_{k}\rangle\langle
M^{(i)}_{l}\rangle\;,$ (18)

which gives the polynomials

$\displaystyle\mathcal{S}_{00}$ $\displaystyle=\mathcal{S}^{2}_{0}-N,$ (19)
$\displaystyle\mathcal{S}_{11}$ $\displaystyle=\mathcal{S}^{2}_{1}-N,$ (20)
$\displaystyle\mathcal{S}_{01}$
$\displaystyle=\mathcal{S}_{0}\mathcal{S}_{1}-(a-b-c+d)$ (21)

Performing the same decomposition as Eq. (II) for three- and four-body
correlators we get

$\displaystyle\mathcal{S}_{000}$
$\displaystyle=\mathcal{S}^{3}_{0}+2\mathcal{S}_{0}-3N\mathcal{S}_{0}$ (22)
$\displaystyle\mathcal{S}_{111}$
$\displaystyle=\mathcal{S}^{3}_{1}+2\mathcal{S}_{1}-3N\mathcal{S}_{1}$ (23)
$\displaystyle\mathcal{S}_{001}$
$\displaystyle=\mathcal{S}_{0}\mathcal{S}_{0}\mathcal{S}_{1}+2\mathcal{S}_{1}-N\mathcal{S}_{1}-2(a-b-c+d)\mathcal{S}_{0}$
(24) $\displaystyle\mathcal{S}_{011}$
$\displaystyle=\mathcal{S}_{0}\mathcal{S}_{1}\mathcal{S}_{1}+2\mathcal{S}_{0}-N\mathcal{S}_{0}-2(a-b-c+d)\mathcal{S}_{1}.$
(25) $\displaystyle\mathcal{S}_{0000}$
$\displaystyle=\mathcal{S}_{0}\mathcal{S}_{0}\mathcal{S}_{0}\mathcal{S}_{0}-6N+3N^{2}-6N\mathcal{S}_{0}\mathcal{S}_{0}+8\mathcal{S}_{0}\mathcal{S}_{0}$
(26) $\displaystyle\mathcal{S}_{1111}$
$\displaystyle=\mathcal{S}_{1}\mathcal{S}_{1}\mathcal{S}_{1}\mathcal{S}_{1}-6N+3N^{2}-6N\mathcal{S}_{1}\mathcal{S}_{1}+8\mathcal{S}_{1}\mathcal{S}_{1}$
(27) $\displaystyle\mathcal{S}_{0001}$
$\displaystyle=\mathcal{S}_{0}\mathcal{S}_{0}\mathcal{S}_{0}\mathcal{S}_{1}-6(a-b-c+d)+3N(a-b-c+d)-3N\mathcal{S}_{0}\mathcal{S}_{1}-3(a-b-c+d)\mathcal{S}_{0}\mathcal{S}_{0}+8\mathcal{S}_{0}\mathcal{S}_{1}$
(28) $\displaystyle\mathcal{S}_{0111}$
$\displaystyle=\mathcal{S}_{0}\mathcal{S}_{1}\mathcal{S}_{1}\mathcal{S}_{1}-6(a-b-c+d)+3N(a-b-c+d)-3N\mathcal{S}_{0}\mathcal{S}_{1}-3(a-b-c+d)\mathcal{S}_{1}\mathcal{S}_{1}+8\mathcal{S}_{0}\mathcal{S}_{1}$
(29) $\displaystyle\mathcal{S}_{0011}$
$\displaystyle=\mathcal{S}_{0}\mathcal{S}_{0}\mathcal{S}_{1}\mathcal{S}_{1}-6N+N^{2}+2(a-b-c+d)^{2}-N\mathcal{S}_{1}\mathcal{S}_{1}-N\mathcal{S}_{0}\mathcal{S}_{0}-4(a-b-c+d)\mathcal{S}_{0}\mathcal{S}_{1}+4\mathcal{S}_{1}\mathcal{S}_{1}+4\mathcal{S}_{0}\mathcal{S}_{0}\;.$
(30)

###  II.1 Proof of the classical bound for $I_{3}$

Here we prove that $I_{3}$ is a valid Bell inequality for arbitrary $N$, by
showing that it cannot be violated by LHVM. The results presented in the
previous section allow us to rewrite the Bell inequality $I_{3}$ of the main
text as a polynomial in the $\\{a,b,c,d\\}$ variables, namely

$\displaystyle I_{3}(a,b,c,d)=8(a-d)(a-d-1)(3c+2d+a-2)+48b(c+d)\;.$ (31)

If $I_{3}(a,b,c,d)\geq 0$ for all integer partitions of $N$ into
$\\{a,b,c,d\\}$, then we must conclude that $I_{3}$ cannot be violated by
LHVM, and thus that it is a valid Bell inequality.

Before starting with the proof, let us note that since $b,c,d\geq 0$, we have

$\displaystyle I_{3}(a,b,c,d)\geq 8(a-d)(a-d-1)(3c+2d+a-2)\;.$ (32)

Let us then begin by proving that

$\displaystyle(a-d)(a-d-1)\geq 0,\;\;\forall(a,d)\in\mathbb{Z}_{\geq
0}^{2}\;.$ (33)

This can be seen as an equation in $e\equiv a-d$ which reads

$\displaystyle e(e-1)=(e-1/2)^{2}-1/4,\;\;e\in\mathbb{Z}\;,$ (34)

which means that $e(e-1)\geq 0$ for all integer $e$: It is strictly positive
if $e>1$ from the sum of squares decomposition and zero if $e=0$ or $e=1$. It
does not matter that the minimum over the reals is $e=-1/2$ for which the
expression evaluates to $-1/4<0$.

Hence, it remains to prove that $3c+2d+a-2\geq 0$ on nonnegative integers.
Note, however, that we do not need to consider the cases where $a-d=0$ or
$a-d-1=0$, since these are prefactors that will make the whole expression
zero.

Since $c\geq 0$, clearly $3c+2d+a-2\geq 2d+a-2$. This last expression can be
negative on the following pairs $(a,d)$:

  * •

$(a,d)=(0,0)$. This is discarded by the premultiplication by $(a-d)$, which
evaluates the whole thing to zero.

  * •

$(a,d)=(1,0)$. This is again discarded by the premultiplication by $(a-d-1)$,
which also eventually evaluates to zero.

  * •

If $a>1$ or $d>0$, then $2d+a-2\geq 0$, so there are no more remaining cases.

All the possibilities are exhausted, thus completing the proof that

$\displaystyle I_{3}(a,b,c,s)\geq 0\quad\forall(a,b,c,d)\in\mathbb{Z}_{\geq
0}^{4}\;.$ (35)

This result implies that $I_{3}$ is a valid Bell inequality, as it shows that
it cannot be violated by LHVM. Following the same approach as the one
presented here, one can prove the validity of any other PIBIs.

##  III II. Writing $\mathcal{S}$ correlators in terms of collective spin
expectation values

Our goal is here to show how PI correlators $\mathcal{S}$ can be measured from
collective spin measurements. First, we parametrise our measurements $M_{k}$
as spin projection measurements, namely

$\displaystyle M^{(i)}_{k}=\vec{u}_{k}\cdot\vec{\sigma}^{(i)}\;$ (36)

with $\vec{u}_{k}\in\mathbb{R}^{3}$ and $|\vec{u_{k}}|^{2}=1$. Here, $k=0,1$
labels one of the two measurement settings considered, and
$\vec{\sigma}=\\{\sigma_{x},\sigma_{y},\sigma_{z}\\}$ is the vector of Pauli
matrices. In the following, we will use the relabeling $\vec{u}_{0}=\vec{n}$
and $\vec{u}_{1}=\vec{m}$.

We consider the case where the same pair of observables are chosen by all
parties, i.e. $M_{k}^{(i)}=M_{k}$ independent of $i$. This allows us to write

$\displaystyle\mathcal{S}_{k}=\sum_{i=1}^{N}\langle
M^{(i)}_{k}\rangle=2\hat{S}_{\vec{u}_{k}}\;,$ (37)

where we have introduced the $N$-particle collective spin operator along
direction $\vec{u}_{k}$ as

$\hat{S}_{\vec{u}_{k}}=\dfrac{1}{2}\sum_{i=1}^{N}\vec{u}_{k}\cdot\vec{\sigma}^{(i)}\;.$
(38)

The one-body correlators thus read

$\displaystyle\mathcal{S}_{0}$
$\displaystyle=2\langle\hat{S}_{\vec{n}}\rangle,$ (39)
$\displaystyle\mathcal{S}_{1}$
$\displaystyle=2\langle\hat{S}_{\vec{m}}\rangle\;.$ (40)

Following the same idea, we can express correlators of any order as collective
spin measurements. For example, for the two-body correlators we can write

$\displaystyle\mathcal{S}_{kl}=\sum_{i,j=1(i\neq j)}^{N}\langle
M_{k}^{(i)}M_{l}^{(j)}\rangle=\sum_{i,j=1}^{N}\langle
M_{k}^{(i)}M_{l}^{(j)}\rangle-\sum_{i=1}^{N}\langle
M_{k}^{(i)}M_{l}^{(i)}\rangle\;.$ (41)

which gives

$\displaystyle\mathcal{S}_{00}$
$\displaystyle=4\langle\hat{S}_{\vec{n}}^{2}\rangle-N,$ (42)
$\displaystyle\mathcal{S}_{11}$
$\displaystyle=4\langle\hat{S}_{\vec{m}}^{2}\rangle-N,$ (43)
$\displaystyle\mathcal{S}_{01}$
$\displaystyle=2\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\rangle-(\vec{m}\cdot\vec{n})N\;.$
(44)

The three-body correlators can be expressed as

$\displaystyle\mathcal{S}_{klm}$ $\displaystyle=\sum^{N}_{i,j,y=1(i\neq j\neq
y)}\langle M_{k}^{(i)}M_{l}^{(j)}M_{m}^{(y)}\rangle$
$\displaystyle=\sum^{N}_{i,j,y=1}\langle
M_{k}^{(i)}M_{l}^{(j)}M_{m}^{(y)}\rangle+2\sum^{N}_{i=1}\langle
M_{k}^{(i)}M_{l}^{(i)}M_{m}^{(i)}\rangle$
$\displaystyle\;-\langle\sum^{N}_{i=1}M^{(i)}_{k}M^{(i)}_{l}\sum^{N}_{y=1}M^{(y)}_{m}\rangle-\langle\sum^{N}_{j=1}M^{(j)}_{l}M^{(j)}_{m}\sum^{N}_{i=1}M^{(i)}_{k}\rangle-\langle\sum^{N}_{y=1}M^{(y)}_{m}M^{(y)}_{k}\sum^{N}_{j=1}M^{(j)}_{l}\rangle\;,$
(45)

which gives

$\displaystyle\mathcal{S}_{000}$
$\displaystyle=8\langle\hat{S}^{3}_{\vec{n}}\rangle+4\langle\hat{S}_{\vec{n}}\rangle-6N\langle\hat{S}_{\vec{n}}\rangle,$
(46) $\displaystyle\mathcal{S}_{111}$
$\displaystyle=8\langle\hat{S}^{3}_{\vec{m}}\rangle+4\langle\hat{S}_{\vec{m}}\rangle-6N\langle\hat{S}_{\vec{m}}\rangle,$
(47) $\displaystyle\mathcal{S}_{001}$
$\displaystyle=\frac{1}{3}\left[8\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\rangle+(4-6N)\langle\hat{S}_{\vec{m}}\rangle+(8-12N)(\vec{n}\cdot\vec{m})\langle\hat{S}_{\vec{n}}\rangle\right],$
(48) $\displaystyle\mathcal{S}_{110}$
$\displaystyle=\frac{1}{3}\left[8\langle\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\rangle+(4-6N)\langle\hat{S}_{\vec{n}}\rangle+(8-12N)(\vec{n}\cdot\vec{m})\langle\hat{S}_{\vec{m}}\rangle\right]\;.$
(49)

Finally, the four-body correlators can be expressed as

$\displaystyle\mathcal{S}_{klmn}$ $\displaystyle=\sum_{i,j,y,h=1(i\neq j\neq
y\neq h)}\langle M^{(i)}_{k}M^{(j)}_{l}M^{(y)}_{m}M^{(h)}_{n}\rangle$
$\displaystyle=\langle\sum^{N}_{i,j,y,h=1}M^{(i)}_{k}M^{(j)}_{l}M^{(y)}_{m}M^{(h)}_{n}\rangle-6\langle\sum^{N}_{i=1}M^{(i)}_{k}M^{(i)}_{l}M^{(i)}_{m}M^{(i)}_{n}\rangle$
$\displaystyle\;+\langle\sum^{N}_{i=1}M^{(i)}_{k}M^{(i)}_{l}\sum^{N}_{y=1}M^{(y)}_{m}M^{(y)}_{n}\rangle+\langle\sum^{N}_{i=1}M^{(i)}_{k}M^{(i)}_{m}\sum^{N}_{j=1}M^{(j)}_{l}M^{(j)}_{n}\rangle+\langle\sum^{N}_{i=1}M^{(i)}_{k}M^{(i)}_{n}\sum^{N}_{j=1}M^{(j)}_{l}M^{(j)}_{m}\rangle$
$\displaystyle\;-\langle\sum^{N}_{i=1}M^{(i)}_{k}M^{(i)}_{l}\sum^{N}_{y=1}M^{(y)}_{m}\sum^{N}_{h=1}M^{(h)}_{n}\rangle-\langle\sum^{N}_{i=1}M^{(i)}_{k}M^{(i)}_{m}\sum^{N}_{j=1}M^{(j)}_{l}\sum^{N}_{h=1}M^{(h)}_{n}\rangle-\langle\sum^{N}_{i=1}M^{(i)}_{k}M^{(i)}_{n}\sum^{N}_{j=1}M^{(j)}_{l}\sum^{N}_{y=1}M^{(y)}_{m}\rangle$
$\displaystyle\;-\langle\sum^{N}_{j=1}M^{(j)}_{l}M^{(j)}_{m}\sum^{N}_{i=1}M^{(i)}_{k}\sum^{N}_{h=1}M^{(h)}_{n}\rangle-\langle\sum^{N}_{j=1}M^{(j)}_{l}M^{(j)}_{n}\sum^{N}_{i=1}M^{(i)}_{k}\sum^{N}_{y=1}M^{(y)}_{m}\rangle-\langle\sum^{N}_{y=1}M^{(y)}_{m}M^{(y)}_{n}\sum^{N}_{i=1}M^{(i)}_{k}\sum^{N}_{j=1}M^{(j)}_{l}\rangle$
$\displaystyle\;+2\langle\sum^{N}_{i=1}M^{(i)}_{k}M^{(i)}_{l}M^{(i)}_{m}\sum^{N}_{h=1}M^{(h)}_{n}\rangle+2\langle\sum^{N}_{j=1}M^{(j)}_{l}M^{(j)}_{m}M^{(j)}_{n}\sum^{N}_{i=1}M^{(i)}_{k}\rangle+2\langle\sum^{N}_{y=1}M^{(y)}_{m}M^{(y)}_{n}M^{(y)}_{k}\sum^{N}_{j=1}M^{(j)}_{l}\rangle+2\langle\sum^{N}_{h=1}M^{(h)}_{n}M^{(h)}_{k}M^{(h)}_{l}\sum^{N}_{y=1}M^{(y)}_{m}\rangle\;.$
(50)

which gives

$\displaystyle\mathcal{S}_{0000}$
$\displaystyle=16\langle\hat{S}^{4}_{\vec{n}}\rangle-6N+3N^{2}-24N\langle\hat{S}^{2}_{\vec{n}}\rangle+32\langle\hat{S}^{2}_{\vec{n}}\rangle,$
(51) $\displaystyle\mathcal{S}_{1111}$
$\displaystyle=16\langle\hat{S}^{4}_{\vec{m}}\rangle-6N+3N^{2}-24N\langle\hat{S}^{2}_{\vec{m}}\rangle+32\langle\hat{S}^{2}_{\vec{m}}\rangle,$
(52) $\displaystyle\mathcal{S}_{0001}$
$\displaystyle=4\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\rangle+(3N^{2}-6N)(\vec{n}\cdot\vec{m})$
$\displaystyle\;+(8-6N)\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\rangle+(16-12N)(\vec{n}\cdot\vec{m})\langle\hat{S}^{2}_{\vec{n}}\rangle,$
(53) $\displaystyle\mathcal{S}_{0111}$
$\displaystyle=4\langle\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\rangle+(3N^{2}-6N)(\vec{n}\cdot\vec{m})$
$\displaystyle\;+(8-6N)\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\rangle+(16-12N)(\vec{n}\cdot\vec{m})\langle\hat{S}^{2}_{\vec{m}}\rangle,$
(54) $\displaystyle\mathcal{S}_{0011}$
$\displaystyle=\frac{1}{6}[16\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\rangle$
$\displaystyle\;-12N+6N^{2}+12N^{2}(\vec{n}\cdot\vec{m})^{2}-24N(\vec{n}\cdot\vec{m})^{2}$
$\displaystyle\;-24N\langle\hat{S}^{2}_{\vec{m}}\rangle-24N\langle\hat{S}^{2}_{\vec{n}}\rangle-48N(\vec{n}\cdot\vec{m})\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\rangle$
$\displaystyle\;+32\langle\hat{S}^{2}_{\vec{m}}\rangle+32\langle\hat{S}^{2}_{\vec{n}}\rangle+64(\vec{n}\cdot\vec{m})\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\rangle]\;.$
(55)

These results allows us to express straightforwardly any PIBI involving
correlators of order $K\leq 4$ in terms of collective spin operator. The
result is a so-called Bell correlation witness, which is an inequality of
practical experimental test. For PIBIs involving correlators of order $K>4$,
expressions similar to the one presented in this section have to be found, a
task which can be significantly simplified by running a computer program.

##  IV III. Expectation values of collective spin high-order moments for OAT
states

For a system of $N$ spin-1/2 particles, the OAT Hamiltonian in terms of the
collective spin operator reads $H=\hbar\chi S_{z}^{2}$. After an adimensional
interaction time $\mu=2\chi t$, a coherent spin state along the $x$-direction
evolves into the state

$\displaystyle|\Phi(\mu)\rangle=\frac{1}{\sqrt{2^{N}}}\sum_{k=0}^{N}\sqrt{\binom{N}{k}}e^{-i(N/2-k)^{2}\mu/2}|k\rangle\;,$
(56)

where $|k\rangle$ denotes the $N$-particle Dicke state with $k$ excitations.

Expectation values of arbitrary moments of the collective spin can be computed
analytically. For a collective spin along direction
$\vec{u}=(u_{x},u_{y},u_{z})$, we find the moments up to order four to be

$\displaystyle\langle\hat{S}_{\vec{u}}\rangle$
$\displaystyle=Su_{x}\cos^{-1+2S}{(\mu/2)},$ (57)
$\displaystyle\langle\hat{S}_{\vec{u}}^{2}\rangle$
$\displaystyle=\frac{1}{4}S\left[(1+2S)(u_{x}^{2}+u_{y}^{2})+2u_{z}^{2}+(-1+2S)\left((u_{x}-u_{y})(u_{x}+u_{y})\cos^{2(-1+S)}{(\mu)}+4u_{y}u_{z}\cos^{-2+2S}{(\mu/2)}\sin{(\mu/2)}\right)\right],$
(58) $\displaystyle\langle\hat{S}_{\vec{u}}^{3}\rangle$
$\displaystyle=\frac{1}{8}S\left[u_{x}\cos^{-3+2S}{(\mu/2)}\left((1-3S+6S^{2})(u_{x}^{2}+u_{y}^{2})-4(2+3(-2+S)S)u_{z}^{2}+2\left((-1+3S)(u_{x}^{2}+u_{y}^{2})+2(1+3(-1+S)S)u_{z}^{2}\right)\cos{(\mu)}\right)\right.$
$\displaystyle\;\left.+(-1+S)(-1+2S)u_{x}(u_{x}^{2}-3u_{y}^{2})\cos^{-3+2S}{(3\mu/2)}+12(-1+S)(-1+2S)u_{x}u_{y}u_{z}\cos^{-3+2S}{(\mu)}\sin{(\mu)}\right],$
(59) $\displaystyle\langle\hat{S}_{\vec{u}}^{4}\rangle$
$\displaystyle=\frac{1}{32}S\left\\{(-1+S+12S^{2}+12S^{3})(u_{x}^{2}+u_{y}^{2})^{2}+8S(-1+6S)(u_{x}^{2}+u_{y}^{2})u_{z}^{2}+8(-1+3S)u_{z}^{4}\right.$
$\displaystyle\;+(-1+2S)\left[(-1+S)(-3+2S)(u_{x}^{4}-6u_{x}^{2}u_{y}^{2}+u_{y}^{4})\cos^{-4+2S}(2\mu)+4(u_{x}-u_{y})(u_{x}+u_{y})\cos^{-4+2S}{(\mu)}\right.$
$\displaystyle\;\times\left(u_{x}^{2}-2Su_{x}^{2}+2S^{2}u_{x}^{2}+u_{y}^{2}-2Su_{y}^{2}+2S^{2}u_{y}^{2}-11u_{z}^{2}+18Su_{z}^{2}-6S^{2}u_{z}^{2}-2u_{x}^{2}\cos{(2\mu)}+3Su_{x}^{2}\cos{(2\mu)}\right.$
$\displaystyle\;\left.-2u_{y}^{2}\cos{(2\mu)}+3Su_{y}^{2}\cos{(2\mu)}+7u_{z}^{2}\cos{(2\mu)}-12Su_{z}^{2}\cos{(2\mu)}+6S^{2}u_{z}^{2}\cos{(2\mu)}\right)$
$\displaystyle\;-8(-1+S)(-3+2S)u_{y}(-3u_{x}^{2}+u_{y}^{2})u_{z}\cos^{-4+2S}{(3\mu/2)}\sin{(3\mu/2)}+8u_{y}u_{z}\cos^{-4+2S}{(\mu/2)}$
$\displaystyle\;\times\left(7u_{x}^{2}\sin{(\mu/2)}-12Su_{x}^{2}\sin{(\mu/2)}+6S^{2}u_{x}^{2}\sin{(\mu/2)}+7u_{y}^{2}\sin{(\mu/2)}-12Su_{y}^{2}\sin{(\mu/2)}+6S^{2}u_{y}^{2}\sin{(\mu/2)}\right.$
$\displaystyle\;-11u_{z}^{2}\sin{(\mu/2)}+18Su_{z}^{2}\sin{(\mu/2)}-6S^{2}u_{z}^{2}\sin{(\mu/2)}-2u_{x}^{2}\sin{(3\mu/2)}+3Su_{x}^{2}\sin{(3\mu/2)}$
$\displaystyle\;\left.\left.\left.-2u_{y}^{2}\sin{(3\mu/2)}+3Su_{y}^{2}\sin{(3\mu/2)}+u_{z}^{2}\sin{(3\mu/2)}-2Su_{z}^{2}\sin{(3\mu/2)}+2S^{2}u_{z}^{2}\sin{(3\mu/2)}\right)\right]\right\\},$
(60)

where $S=N/2$ is total spin.

These expressions allow us to conveniently compute the violation of Bell
correlation witnesses for OAT states, as a function of the chosen measurement
directions.

##  V IV. Writing collective spin expectation values in terms of $\mathcal{S}$
correlators

To construct experimentally practical Bell operators, our approach consists of
starting from a set of convenient measurements, and derive from them the
structure of the associated Bell inequalities. This is enabled by the fact
that expectation values of a collective spin operator can be written in terms
of the permutationally invariant $\mathcal{S}$ correlators. To see this, note
first that the results presented in Supplementary Sec.II can be inverted to
write

$\displaystyle\langle\hat{S}_{\vec{n}}\rangle=\frac{\mathcal{S}_{0}}{2},$ (61)
$\displaystyle\langle\hat{S}_{\vec{m}}\rangle=\frac{\mathcal{S}_{1}}{2},$ (62)
$\displaystyle\langle\hat{S}^{2}_{\vec{n}}\rangle=\frac{\mathcal{S}_{00}+N}{4},$
(63)
$\displaystyle\langle\hat{S}^{2}_{\vec{m}}\rangle=\frac{\mathcal{S}_{11}+N}{4},$
(64)
$\displaystyle\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\rangle=\frac{\mathcal{S}_{01}+(\vec{n}\cdot\vec{m})N}{2},$
(65)
$\displaystyle\langle\hat{S}^{3}_{\vec{n}}\rangle=\frac{\mathcal{S}_{000}+(3N-2)\mathcal{S}_{0}}{8},$
(66)
$\displaystyle\langle\hat{S}^{3}_{\vec{m}}\rangle=\frac{\mathcal{S}_{111}+(3N-2)\mathcal{S}_{1}}{8},$
(67)
$\displaystyle\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\rangle=\frac{3\mathcal{S}_{001}-(2-3N)\mathcal{S}_{1}-(4-6N)(\vec{n}\cdot\vec{m})\mathcal{S}_{0}}{8},$
(68)
$\displaystyle\langle\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\rangle=\frac{3\mathcal{S}_{011}-(2-3N)\mathcal{S}_{0}-(4-6N)(\vec{n}\cdot\vec{m})\mathcal{S}_{1}}{8},$
(69)
$\displaystyle\langle\hat{S}^{4}_{\vec{n}}\rangle=\frac{\mathcal{S}_{0000}+2(3N-4)(\mathcal{S}_{00}+N)-3(N-2)N}{16},$
(70)
$\displaystyle\langle\hat{S}^{4}_{\vec{m}}\rangle=\frac{\mathcal{S}_{1111}+2(3N-4)(\mathcal{S}_{11}+N)-3(N-2)N}{16},$
(71)
$\displaystyle\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\rangle$
$\displaystyle\qquad=\frac{\mathcal{S}_{0001}-3(N-2)N(\vec{n}\cdot\vec{m})+(3N-4)(\vec{n}\cdot\vec{m})(\mathcal{S}_{00}+N)-(4-3N)(\mathcal{S}_{01}+N(\vec{n}\cdot\vec{m}))}{4},$
(72)
$\displaystyle\langle\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\rangle$
$\displaystyle\qquad=\frac{\mathcal{S}_{0111}-3(N-2)N(\vec{n}\cdot\vec{m})+(3N-4)(\vec{n}\cdot\vec{m})(\mathcal{S}_{11}+N)-(4-3N)(\mathcal{S}_{01}+N(\vec{n}\cdot\vec{m}))}{4},$
(73)
$\displaystyle\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\rangle$
$\displaystyle\qquad=\frac{3\mathcal{S}_{0011}+(3N-4)(\mathcal{S}_{00}+\mathcal{S}_{11}+4(\vec{n}\cdot\vec{m})\mathcal{S}_{01})+N(3N-2)\left(2(\vec{n}\cdot\vec{m})^{2}+1\right)}{8}\;.$
(74)

Using the (linearly independent) measurement direction $\vec{m}$ and
$\vec{n}$, we can write any measurement direction $\vec{a}$ lying on the
$\\{\vec{m},\vec{n}\\}$ plane as

$\displaystyle\vec{a}=\alpha\vec{m}+\beta\vec{n}$ (75)

with normalized weights $|\alpha|^{2}+|\beta|^{2}=1$, and
$\alpha,\beta\in\mathbb{R}^{3}$.

Expectation values for the first four moments of the collective spin operator
along some direction $\vec{a}$ are thus

$\displaystyle\langle\hat{S}_{\vec{a}}\rangle$
$\displaystyle=\alpha\langle\hat{S}_{\vec{m}}\rangle+\beta\langle\hat{S}_{\vec{n}}\rangle$
$\displaystyle=\frac{\beta}{2}\mathcal{S}_{0}+\frac{\alpha}{2}\mathcal{S}_{1},$
(76) $\displaystyle\langle\hat{S}^{2}_{\vec{a}}\rangle$
$\displaystyle=\alpha^{2}\langle\hat{S}^{2}_{\vec{m}}\rangle+\beta^{2}\langle\hat{S}^{2}_{\vec{n}}\rangle+\alpha\beta\langle\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\rangle$
$\displaystyle=\frac{\beta^{2}}{4}\mathcal{S}_{00}+\frac{\alpha^{2}}{4}\mathcal{S}_{11}+\frac{\alpha\beta}{2}\mathcal{S}_{01}+\left(\alpha^{2}+\beta^{2}+2\alpha\beta(\vec{n}\cdot\vec{m})\right)\frac{N}{4},$
(77) $\displaystyle\langle\hat{S}^{3}_{\vec{a}}\rangle$
$\displaystyle=\alpha^{3}\langle\hat{S}^{3}_{\vec{m}}\rangle+\beta^{3}\langle\hat{S}^{3}_{\vec{n}}\rangle+\alpha^{2}\beta\langle\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\rangle+\alpha\beta^{2}\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\rangle$
$\displaystyle=\frac{\beta^{3}}{8}\mathcal{S}_{000}+\frac{\alpha^{3}}{8}\mathcal{S}_{111}+\frac{3\alpha\beta^{2}}{8}\mathcal{S}_{001}+\frac{3\alpha^{2}\beta}{8}\mathcal{S}_{011}$
$\displaystyle\;+\frac{(3N-2)}{8}\left(\beta^{3}+\alpha^{2}\beta+2\alpha\beta^{2}(\vec{n}\cdot\vec{m})\right)\mathcal{S}_{0}+\frac{(3N-2)}{8}\left(\alpha^{3}+\alpha\beta^{2}+2\alpha^{2}\beta(\vec{n}\cdot\vec{m})\right)\mathcal{S}_{1}$
(78) $\displaystyle\langle\hat{S}_{\vec{a}}^{4}\rangle$
$\displaystyle=\alpha^{4}\langle\hat{S}^{4}_{\vec{m}}\rangle+\beta^{4}\langle\hat{S}^{4}_{\vec{n}}\rangle+\alpha^{3}\beta\langle\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\rangle$
$\displaystyle\;+\alpha\beta^{3}\langle\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\rangle$
$\displaystyle\;+\alpha^{2}\beta^{2}\langle\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}+\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}+\hat{S}_{\vec{n}}\hat{S}_{\vec{m}}\hat{S}_{\vec{m}}\hat{S}_{\vec{n}}\rangle$
$\displaystyle=\frac{\beta^{4}}{16}\mathcal{S}_{0000}+\frac{\alpha^{4}}{16}\mathcal{S}_{1111}+\frac{\alpha\beta^{3}}{4}\mathcal{S}_{0001}+\frac{\alpha^{3}\beta}{4}\mathcal{S}_{0111}+\frac{3\alpha^{2}\beta^{2}}{8}\mathcal{S}_{0011}$
$\displaystyle\;+(3N-4)\left(\frac{\beta^{4}}{8}+\frac{\alpha\beta^{3}}{4}(\vec{n}\cdot\vec{m})+\frac{\alpha^{2}\beta^{2}}{8}\right)\mathcal{S}_{00}$
$\displaystyle\;+(3N-4)\left(\frac{\alpha^{4}}{8}+\frac{\alpha^{3}\beta}{4}(\vec{n}\cdot\vec{m})+\frac{\alpha^{2}\beta^{2}}{8}\right)\mathcal{S}_{11}$
$\displaystyle\;+(3N-4)\left(\frac{\alpha^{3}\beta}{4}+\frac{\alpha\beta^{3}}{4}+\frac{\alpha^{2}\beta^{2}}{2}(\vec{n}\cdot\vec{m})\right)\mathcal{S}_{01}$
$\displaystyle\;+N(3N-2)\left(\frac{\alpha^{4}}{16}+\frac{\beta^{4}}{16}+\frac{\alpha^{3}\beta}{4}(\vec{n}\cdot\vec{m})+\frac{\alpha\beta^{3}}{4}(\vec{n}\cdot\vec{m})+\frac{\alpha^{2}\beta^{2}}{8}(2(\vec{n}\cdot\vec{m})^{2}+1)\right).$
(79)

These results allow us to derive the form of the PIBI associated to a specific
Bell operator.

##  VI V. Families of PIBIs involving third-order correlators

Besides the PIBI $I_{3}$ presented in the main text, we also find the
following PIBIs involving third-order correlators. See Fig. 4 for a plot of
their relative quantum violation as a function of $N$.

$\displaystyle I_{3}^{(1)}$
$\displaystyle\equiv-24(N-1)\mathcal{S}_{0}+3(N+2)\mathcal{S}_{00}-6(N-2)\mathcal{S}_{01}+3(N-2)\mathcal{S}_{11}-2\mathcal{S}_{000}+3\mathcal{S}_{001}-\mathcal{S}_{111}+12N(N-1)\geq
0,$ $\displaystyle I_{3}^{(2)}$
$\displaystyle\equiv-4(3N-1)\mathcal{S}_{0}-4(N-1)\mathcal{S}_{1}+(N+6)\mathcal{S}_{00}+2N\mathcal{S}_{01}+(N-2)\mathcal{S}_{11}-\mathcal{S}_{000}-2\mathcal{S}_{001}-\mathcal{S}_{011}+4N(N+1)\geq
0,$ $\displaystyle I_{3}^{(3)}$
$\displaystyle\equiv-8N\mathcal{S}_{0}-8(N-1)\mathcal{S}_{1}+(N+2)\mathcal{S}_{00}+2(N+2)\mathcal{S}_{01}+(N-2)\mathcal{S}_{11}-\mathcal{S}_{000}-2\mathcal{S}_{001}-\mathcal{S}_{011}+4N(N+1)\geq
0,$ $\displaystyle I_{3}^{(4)}$
$\displaystyle\equiv-8(N-1)\mathcal{S}_{0}-8\mathcal{S}_{1}+(N-2)\mathcal{S}_{00}-2(N-6)\mathcal{S}_{01}+(N-6)\mathcal{S}_{11}-\mathcal{S}_{001}+2\mathcal{S}_{011}-\mathcal{S}_{111}+4N(N+1)\geq
0,$ $\displaystyle I_{3}^{(5)}$
$\displaystyle\equiv-8(N-1)\mathcal{S}_{0}+(N+2)\mathcal{S}_{00}-2(N-2)\mathcal{S}_{01}+(N-2)\mathcal{S}_{11}-\mathcal{S}_{000}+2\mathcal{S}_{001}-\mathcal{S}_{011}+4N(N-1)\geq
0,$ $\displaystyle I_{3}^{(6)}$
$\displaystyle\equiv-4(N-1)\mathcal{S}_{0}-4(N-1)\mathcal{S}_{1}+(N-2)\mathcal{S}_{00}-2(N-4)\mathcal{S}_{01}+(N-2)\mathcal{S}_{11}-\mathcal{S}_{000}+2\mathcal{S}_{001}-\mathcal{S}_{011}+4N(N-1)\geq
0,$ $\displaystyle I_{3}^{(7)}$
$\displaystyle\equiv-4(N-1)\mathcal{S}_{0}-4(N-1)\mathcal{S}_{1}+(N-2)\mathcal{S}_{00}+2N\mathcal{S}_{01}+(N-2)\mathcal{S}_{11}-\mathcal{S}_{000}-2\mathcal{S}_{001}-\mathcal{S}_{011}+4N(N-1)\geq
0,$ $\displaystyle I_{3}^{(8)}$
$\displaystyle\equiv(N-6)\mathcal{S}_{00}-2(N-2)\mathcal{S}_{01}+(N-2)\mathcal{S}_{11}-\mathcal{S}_{000}+2\mathcal{S}_{001}-\mathcal{S}_{011}+4N(N-1)\geq
0,$ $\displaystyle I_{3}^{(9)}$ $\displaystyle\equiv
3(N-4)\mathcal{S}_{00}-6N\mathcal{S}_{01}+3N\mathcal{S}_{11}-2\mathcal{S}_{000}+3\mathcal{S}_{001}-\mathcal{S}_{111}+12N(N-1)\geq
0,$ $\displaystyle I_{3}^{(10)}$
$\displaystyle\equiv-12(N-1)\mathcal{S}_{0}-4(N-7)\mathcal{S}_{1}+(N-2)\mathcal{S}_{00}+2(N-8)\mathcal{S}_{01}+(N-10)\mathcal{S}_{11}+\mathcal{S}_{001}+2\mathcal{S}_{011}+\mathcal{S}_{111}+4N(N+5)\geq
0,$ $\displaystyle I_{3}^{(11)}$
$\displaystyle\equiv-12(N-1)\mathcal{S}_{0}+12(N+1)\mathcal{S}_{1}+(N-2)\mathcal{S}_{00}-2(N+4)\mathcal{S}_{01}+(N+6)\mathcal{S}_{11}+\mathcal{S}_{001}-2\mathcal{S}_{011}+\mathcal{S}_{111}+4N(N+5)\geq
0,$ $\displaystyle I_{3}^{(12)}$ $\displaystyle\equiv
8(N-1)\mathcal{S}_{0}+8(2N+1)\mathcal{S}_{1}+(N-2)\mathcal{S}_{00}+2(N+2)\mathcal{S}_{01}+(N+10)\mathcal{S}_{11}+\mathcal{S}_{001}+2\mathcal{S}_{011}+\mathcal{S}_{111}+4N(N+5)\geq
0,$ $\displaystyle I_{3}^{(13)}$
$\displaystyle\equiv-16(N-1)\mathcal{S}_{0}-8(N-7)\mathcal{S}_{1}+(N-2)\mathcal{S}_{00}+2(N-10)\mathcal{S}_{01}+(N-14)\mathcal{S}_{11}+\mathcal{S}_{001}+2\mathcal{S}_{011}+\mathcal{S}_{111}+4N(N+11)\geq
0,$ $\displaystyle I_{3}^{(14)}$
$\displaystyle\equiv-16(N-1)\mathcal{S}_{0}+16(N+2)\mathcal{S}_{1}+(N-2)\mathcal{S}_{00}-2(N+6)\mathcal{S}_{01}+(N+10)\mathcal{S}_{11}+\mathcal{S}_{001}-2\mathcal{S}_{011}+\mathcal{S}_{111}+4N(N+11)\geq
0,$ $\displaystyle I_{3}^{(15)}$ $\displaystyle\equiv
12(N-1)\mathcal{S}_{0}+4(5N+7)\mathcal{S}_{1}+(N-2)\mathcal{S}_{00}+2(N+4)\mathcal{S}_{01}+(N+14)\mathcal{S}_{11}+\mathcal{S}_{001}+2\mathcal{S}_{011}+\mathcal{S}_{111}+4N(N+11)\geq
0,$ $\displaystyle I_{3}^{(16)}$
$\displaystyle\equiv(20N-44)\mathcal{S}_{0}+(2N+2)\mathcal{S}_{00}+(-4N+16)\mathcal{S}_{01}+(2N-8)\mathcal{S}_{11}+\mathcal{S}_{000}-\mathcal{S}_{001}-\mathcal{S}_{011}+\mathcal{S}_{111}+(10N^{2}-34N+48)\geq
0,$ $\displaystyle I_{3}^{(17)}$ $\displaystyle\equiv
24(N-3)\mathcal{S}_{0}+24(N-3)\mathcal{S}_{1}+5(N-2)\mathcal{S}_{00}-2(5N-22)\mathcal{S}_{01}+5(N-2)\mathcal{S}_{11}+2\mathcal{S}_{000}-\mathcal{S}_{001}-4\mathcal{S}_{011}+3\mathcal{S}_{111}$
$\displaystyle\;+24(N^{2}-5N+8)\geq 0.$ (80)

Figure 4: Maximum relative quantum violation of the PIBIs $I_{3}^{(i)}$,
compared to $I_{3}$ and $I_{2}$ of our main text.

##  VII VI. PIBI involving fourth-order correlators

If we only take into account the second and the fourth order correlators, by
finding the corresponding local polytope we find the PIBI

$\displaystyle
I_{4}=24(N-1)\mathcal{S}_{00}+48(N-1)\mathcal{S}_{01}+24(N-3)\mathcal{S}_{11}+\mathcal{S}_{0000}+4\mathcal{S}_{0001}+6\mathcal{S}_{0011}+4\mathcal{S}_{0111}+\mathcal{S}_{1111}+48N(N-1)\geq
0.$ (81)

In Fig. 5 we compare its maximum relative quantum violation to the one of
$I_{2}$ and $I_{3}$ of the main text, as a function of $N$.

Figure 5: Maximum relative quantum violation of the fourth-order PIBI $I_{4}$,
compared to $I_{3}$ and $I_{2}$ of our main text.

According to Eq. (79), when $\alpha=\beta=\frac{1}{\sqrt{2}}$ we have
$\vec{a}=(\vec{n}+\vec{m})/\sqrt{2}$ and

$\displaystyle\mathcal{S}_{0000}+4\mathcal{S}_{0001}+6\mathcal{S}_{0011}+4\mathcal{S}_{0111}+\mathcal{S}_{1111}$
$\displaystyle=64\langle\hat{S}^{4}_{\vec{a}}\rangle-4(3N-4)(1+(\vec{n}\cdot\vec{m}))(\mathcal{S}_{00}+\mathcal{S}_{11}+2\mathcal{S}_{01})-N(3N-2)\left[2+8(\vec{n}\cdot\vec{m})+2\left(2(\vec{n}\cdot\vec{m})^{2}+1\right)\right]$
(82)

From this, we can rewrite $I_{4}$ as

$\displaystyle I_{4}$
$\displaystyle=64\langle\hat{S}^{4}_{\vec{a}}\rangle+\left[-4(3N-4)(\vec{n}\cdot\vec{m})+4(3N-2)\right]\mathcal{S}_{00}+\left[-4(3N-4)(\vec{n}\cdot\vec{m})+4(3N-14)\right]\mathcal{S}_{11}$
$\displaystyle\;+\left[-8(3N-4)(\vec{n}\cdot\vec{m})+8(3N-2)\right]\mathcal{S}_{01}+\left[-N(3N-2)\left[2+8(\vec{n}\cdot\vec{m})+2\left(2(\vec{n}\cdot\vec{m})^{2}+1\right)\right]+48N(N-1)\right]$
$\displaystyle=64\langle\hat{S}^{4}_{\vec{a}}\rangle-192\langle\hat{S}^{2}_{\vec{m}}\rangle+32[-\sqrt{2}(3N-4)(\vec{a}\cdot\vec{m})+6(N-1)]\langle\hat{S}^{2}_{\vec{a}}\rangle+12(\sqrt{2}(\vec{a}\cdot\vec{m})-2)\left[\sqrt{2}(\vec{a}\cdot\vec{m})(N-2)-2N\right]N$
(83)

which shows that only two measurement directions, namely $\vec{m}$ and
$\vec{a}$, are needed to test $I_{4}$ experimentally.

##  VIII VII. SDP method

In an $N$-party scenario with two measurement settings per party and two
possible outputs per measurement, there exist four local deterministic
strategies (LDSs) per party. For the characterization of permutationally
invariant inequalities, the vertices of the local polytope are naturally
projected onto the symmetric space via permutation of the correlators. Such a
projection onto the symmetric space induces a natural classification of the
local deterministic strategies in terms of the partitions of $N$ in four
elements, which we denote $\\{a,b,c,d\\}$, corresponding to Eq. (II). These
represent how many parties corresponding to each strategy. Hence, $a$, $b$,
$c$, and $d$ are non-negative integers and the sum of them is $a+b+c+d=N$. $a$
can be thought of as the amount of parties that output the first outcome on
the first measurement setting and the first outcome on the second measurement
setting. $b$ the number of parties that output the first and second outcomes
on the first and second measurement settings, respectively, and so on.

In the case where at most third-order correlators are involved, the correlator
ensemble is
$\vec{\mathcal{S}}_{3}\equiv(\mathcal{S}_{0},\mathcal{S}_{1},\mathcal{S}_{00},\mathcal{S}_{01},\mathcal{S}_{1},\mathcal{S}_{000},\mathcal{S}_{001},\mathcal{S}_{011},\mathcal{S}_{111})$.
Furthermore, since correlation functions factorize at each LDS [37], this
factorization induces another factorization at the level of symmetric
correlators, when evaluated at projected vertices [12, 38]. For example, it is
immediate to see that ${\cal S}_{00}={\cal S}_{0}^{2}-N$ or that ${\cal
S}_{01}={\cal S}_{0}{\cal S}_{1}-{\cal Z}$, where ${\cal Z}=a-b-c+d$. This
allows us to express any correlator in terms of linear polynomials in
$a,b,c,d$. More precisely, let us consider the following transformation

$\left(\begin{array}[]{c}N\\\ {\cal S}_{1}\\\ {\cal S}_{0}\\\ {\cal
Z}\end{array}\right)=2H^{\otimes 2}\left(\begin{array}[]{c}a\\\ b\\\ c\\\
d\end{array}\right)=\left(\begin{array}[]{cccc}1&1&1&1\\\ 1&-1&1&-1\\\
1&1&-1&-1\\\ 1&-1&-1&1\end{array}\right)\left(\begin{array}[]{c}a\\\ b\\\ c\\\
d\end{array}\right),$ (84)

where $H$ is the Hadamard gate. However, the description in terms of the
explicit parameters $a,b,c,d$ is not totally convenient for our purposes, nor
is the one in terms of $N,{\cal S}_{1},{\cal S}_{0},{\cal Z}$. The projected
symmetric polytope is embedded in a nine-dimensional affine space and, if
$a,b,c,d$ were real, continuous parameters, these would form a four-
dimensional manifold in that space. The natural coordinates for the projected
symmetric polytope are therefore the symmetric correlators, and $N$ and ${\cal
Z}$ are none of them, but just artifacts of the parameterization. However, we
can relate them to coordinate axis in the polytope affine space, since
$N={\cal S}_{0}^{2}-{\cal S}_{00}$ and ${\cal Z}={\cal S}_{0}{\cal
S}_{1}-{\cal S}_{01}$.

By inverting the above relations, we can transform the inequality constraints
$a,b,c,d\geq 0$ into the following respective polynomial inequalities
$g_{1},\ldots g_{4}$ in the variables of $\vec{\mathcal{S}}_{3}$:

$\displaystyle\left(\begin{matrix}g_{1}\\\ g_{2}\\\ g_{3}\\\
g_{4}\end{matrix}\right)=\frac{1}{4}\left(\begin{matrix}(\mathcal{S}_{0}^{2}-\mathcal{S}_{00})+\mathcal{S}_{0}+\mathcal{S}_{1}+(\mathcal{S}_{0}\mathcal{S}_{1}-\mathcal{S}_{01})\\\
(\mathcal{S}_{0}^{2}-\mathcal{S}_{00})+\mathcal{S}_{0}-\mathcal{S}_{1}-(\mathcal{S}_{0}\mathcal{S}_{1}-\mathcal{S}_{01})\\\
(\mathcal{S}_{0}^{2}-\mathcal{S}_{00})-\mathcal{S}_{0}+\mathcal{S}_{1}-(\mathcal{S}_{0}\mathcal{S}_{1}-\mathcal{S}_{01})\\\
(\mathcal{S}_{0}^{2}-\mathcal{S}_{00})-\mathcal{S}_{0}-\mathcal{S}_{1}+(\mathcal{S}_{0}\mathcal{S}_{1}-\mathcal{S}_{01})\end{matrix}\right)\geq
0,$ (85)

where $g_{i}$ is shorthand of $g_{i}(\vec{\mathcal{S}}_{3})$.

At this stage, we have the tools to take the next step for approximating the
local polytope with a continuous manifold. We first consider
$\\{a,b,c,d\\}\in\mathbb{R}_{\geq 0}$ instead of
$\\{a,b,c,d\\}\in\mathbb{Z}_{\geq 0}$, and relabel them as non-negative
constraints in terms of correlators $g_{i}(\vec{\mathcal{S}}_{3})\geq 0$.
These give a so-called semialgebraic set in the $\vec{\mathcal{S}}_{3}$
coordinates. Semialgebraic sets are defined by polynomial equality and
inequality constraints and they are central objects in convex optimization
[34]. For our purposes, the constraints are the polynomial inequalities
$g_{i}(\vec{\mathcal{S}}_{3})\geq 0$, together with the polynomial equalities
such as ${\cal S}_{000}={\cal
S}_{0}^{3}-(3(\mathcal{S}_{0}^{2}-\mathcal{S}_{00})-2){\cal S}_{0}$,
${\mathcal{S}}_{001}=\mathcal{S}_{0}\mathcal{S}_{0}\mathcal{S}_{1}+2\mathcal{S}_{1}-(\mathcal{S}_{0}^{2}-\mathcal{S}_{00})\mathcal{S}_{1}-2(\mathcal{S}_{0}\mathcal{S}_{1}-\mathcal{S}_{01})\mathcal{S}_{0}$,
etc. (cf. Eqs. (22)-(25)).

The characterization of the convex hull of a discrete set of points
$\\{\vec{\mathcal{S}}_{3}:a,b,c,d\in{\mathbb{Z}},a,b,c,d\geq 0,a+b+c+d=N\\}$
is now relaxed to the characterization of the convex hull of the semialgebraic
set $\\{\vec{\mathcal{S}}_{3}:a,b,c,d\in{\mathbb{R}},a,b,c,d\geq
0,a+b+c+d=N\\}$, which we denote $\overline{CH}({\mathbb{S}})$. By
construction, the latter contains the former. The membership problem in the
convex hull of a semialgebraic set is a subject of intense research in
optimization [34, 35, 36]. Since convex hulls are intersections of half-
spaces, determined by linear inequalities, the hardness of the problem lies in
proving that a linear polynomial is non-negative when evaluated on a
semialgebraic set. However, the question in its full generality is too
complex, so a second relaxation is needed.

The second relaxation is to express the polynomial as a sum of squares so that
the membership problem can be efficiently solved via SDP method [35, 36].
However, the sum of squares has to take into account the polynomial equality
and inequality constraints.

For the equality constraints we introduce the ideal $I$ generated by a set of
polynomials $f_{i}$

$\displaystyle
I=\left\\{\sum_{i}f_{i}p_{i}\;\text{s.t.}\;p_{i}\in\mathbb{R}[\vec{\mathcal{S}}_{3}]\right\\}\subseteq\mathbb{R}[\vec{\mathcal{S}}_{3}],$
(86)

where $f_{i}=f_{i}(\vec{\mathcal{S}}_{3})=0$ is a set of equations associated
to the correlators (such as ${\cal S}_{000}-{\cal
S}_{0}^{3}+(3(\mathcal{S}_{0}^{2}-\mathcal{S}_{00})-2){\cal S}_{0}$). Since
the $f_{i}$ evaluate to zero in the semialgebraic set, every element in $I$
also evaluates to zero. This allows us to tremendously increase expressivity
of the linear polynomials we want to show define half-spaces containing the
semialgebraic set. In particular, $f_{i}$ allow to change the degree of such
polynomials, allowing for sum-of-squares polynomials to be equivalent to e.g.
linear ones when evaluated on that set.

For the inequality constraints, we proceed in a similar way: we can use any of
the $g_{i}$ as non-negative elements in our certificate construction beacause
they are, by construction, non-negative on the semialgebraic set to
characterize.

All in all, this allows us to prove that a linear polynomial
$l(\vec{\mathcal{S}}_{3})\geq 0$ on $\overline{CH}({\mathbb{S}})$ via the
following strategy: If we can express

$l(\vec{\mathcal{S}}_{3})=\sum_{i}g_{i}(\vec{\mathcal{S}}_{3})\sigma_{i}(\vec{\mathcal{S}}_{3})+\sigma_{0}(\vec{\mathcal{S}}_{3})\mod
I,$ (87)

where $\sigma_{i}(\vec{\mathcal{S}}_{3})$ are sum-of-squares polynomials
modulo the ideal $I$ for $0\leq i\leq 4$ (i.e., there exits a set of
polynomials $p_{j}$ such that $\sigma_{i}+\sum_{j}p_{j}f_{j}$ is a sum-of-
squares), then this is a syntactic proof that every element in
$\overline{CH}({\mathbb{S}})$ satisfies $l\geq 0$.

The condition in Eq. (87) can be efficiently computed, if the degree of
$p_{j}$ is bounded, yielding a hierarchy of relaxations. Clearly, the higher
the degree allowed for $p_{j}$, the more possibilities one can express a
particular $l$ as a sum-of-squares modulo $I$, but also the higher the
computational cost. The membership problem in $\overline{CH}({\mathbb{S}})$
can be then relaxed through a suitable modification Lasserre’s moment method
[39, 35, 38].

We focus on the very first level of that hierarchy (both on $\sigma_{0}$ and
$\sigma_{i}$). To this end, we construct a $10\times 10$ moment matrix
$\Gamma_{i}=g_{i}\vec{b}\cdot\vec{b}^{T}$ mod $I$, with the vector
$\vec{b}=(1,\mathcal{S}_{0},\cdots,\mathcal{S}_{111})^{T}$ and linearize its
components. Since we have four inequality constraints also to treat at the
first level, this yields a $50\times 50$ block-diagonal matrix
$\tilde{\Gamma}=\displaystyle\bigoplus_{i=0}^{4}\Gamma_{i}$. Now, some entries
of $\tilde{\Gamma}$ correspond to the correlators that are measurable in the
experiment, yielding a data point denoted $\vec{\mathcal{S}}_{3}^{*}$. If the
semidefinite program

$\begin{array}[]{cccc}\min_{\tilde{\Gamma}}&0&\\\
\mathrm{s.t.}&\tilde{\Gamma}&\succeq&0\\\ &\tilde{\Gamma}_{00}&=&1\\\
&\tilde{\Gamma}_{0i}&=&(\vec{\mathcal{S}}_{3}^{*})_{i}\\\
&\tilde{\Gamma}_{ij}&=&h(\tilde{\Gamma}),\end{array}$ (88)

where $h$ encodes the constraints imposed by the ideal $I$. If the SDP (88) is
feasible, then $\vec{\mathcal{S}}_{3}^{*}$ belongs to the spectrahedron
defined by the feasible set of SDP (88). In that case, the SDP does not
decide, the correlations $\vec{\mathcal{S}}_{3}^{*}$ might admit a LHVM or be
nonlocal. However, if the SDP (88) is not feasible, its dual yields a
certificate of infeasibility. Geometrically speaking, that is a hyperplane
separating the outer approximation of $\overline{CH}({\mathbb{S}})$ given by
the spectrahedron defined by the feasible set of the SDP (88) and the exterior
point $\vec{\mathcal{S}}_{3}^{*}$. Since the spectrahedron contains
$\overline{CH}({\mathbb{S}})$, infeasibility of SDP (88) implies that
$\vec{\mathcal{S}}_{3}^{*}$ is nonlocal. Interestingly, the dual variables
associated to the equality constraints
$\tilde{\Gamma}_{0i}=(\vec{\mathcal{S}}_{3}^{*})_{i}$ yield in that case the
associated coefficients of the symmetric 3-body Bell inequality that detects
$(\vec{\mathcal{S}}_{3}^{*})_{i}$, the dual variable associated to the
constraint $\tilde{\Gamma}_{00}=1$ its classical bound and the dual variable
associated to $\tilde{\Gamma}\succeq 0$ the sum of squares modulo $I$ of the
form (87).

##  IX VIII. Quantitative analysis of high-order PIBI violation and non-
Gaussian properties

To quantitatively investigate the advantage of high-order PIBIs on revealing
Bell nonlocality of non-Gaussian states, in this section we investigate the
relation between PIBIs violation and some measure of spin state non-Gaussian
properties.

First, we look at the excess kurtosis, which is commonly used as a
quantification for non-Gaussianity in continuous-variable systems [40, 41].
For spin states, this quantity can be obtained by optimizing the measurement
direction $\vec{a}=(\phi,\theta)$ of the collective spin $\hat{S}_{\vec{a}}$,
namely

$\displaystyle
K_{\text{ex}}=\min_{\phi,\theta}\left(\frac{\langle\left(\hat{S}_{\vec{a}}-\langle\hat{S}_{\vec{a}}\rangle\right)^{4}\rangle}{\langle\left(\hat{S}_{\vec{a}}-\langle\hat{S}_{\vec{a}}\rangle\right)^{2}\rangle^{2}}-3\right)\;,$
(89)

where $K_{\text{ex}}<0$ indicates non-Gaussianity of the state.

In addition, there are classes of non-Gaussian quantum states that are
associated to Wigner functions with negative regions. For this reason, we also
look at the state’s Wigner negativity, which is defined as the double volume
of the negative part of the Wigner function. For spin states this is [42]

$\mathcal{N}=\frac{1}{2}\left(\frac{2j+1}{4\pi}\int_{\theta=0}^{\pi}\int_{\phi=0}^{2\pi}|W_{\rho}(\theta,\phi)|\sin\theta
d\theta d\phi-1\right)\;,$ (90)

where $W_{\rho}(\theta,\phi)$ is the Wigner function value at point
$(\theta,\phi)$ on the generalized Bloch sphere, and $j=N/2$ is the total
spin. $\mathcal{N}>0$ is an indicator of non-classicality [43], and it has
also been proved to be a necessary ingredient in quantum computation to
outperform classical devices [44].

In Figure 6 we show for OAT states a comparison between the range over which
$I_{2}$ and $I_{3}$ are violated, and the level of non-Gaussian properties of
the states. Fig. 6(a) is a zoom of Fig. 2 in the main text, where we compare
relative quantum violation between PIBI $I_{2}$ and $I_{3}$ with the increase
of OAT interaction $\mu$. The dashed vertical lines indicate the critical
values of $\mu$ where the PIBI violation vanishes. Non-Gaussian properties of
OAT states as a function of $\mu$ are illustrated by excess kurtosis in Fig.
6(b) and Wigner negativity in Fig. 6(c), where it is evident an increase of
non-Gaussianity and non-classicality with $\mu$. By looking at the vertical
dashed lines it is evident that $I_{3}$ outperforms $I_{2}$ in nonlocality
detection, and especially that $I_{3}$ allows us to detect nonlocality in
states with higher levels of non-Gaussianity. This is expected from the fact
that increasing the order of the correlators involved in the PIBI allows for
detecting correlations that do not manifest in lower order correlators.

Figure 6: Comparison between $I_{2}$ and $I_{3}$ violation and non-Gaussian
properties of OAT states with $N=50$. As a function of the OAT interaction
time $\mu$ we present: (a) relative quantum violation
$Q_{V}^{N}/\beta_{V}^{N}$ of the 2nd-order PIBI $I_{2}$ (blue) and 3rd-order
PIBI $I_{3}$ (orange), (b) excess kurtosis $K_{ex}$, and (c) Wigner
negativity. Dashed vertical lines indicate the critical values of $\mu$ where
the PIBIs violation vanishes.

Furthermore, for the highly non-Gaussian state maximising the quantum
violation of $\hat{I}_{4}$ (whose Wigner function is presented in Fig. 3 in
the main text), we find an excess kurtosis of $K_{ex}=-1.94$ and a Wigner
negativity of $\mathcal{N}=1.41$. Note that these numbers are much larger than
the one for the OAT states detected by $I_{2}$ (see Fig. (6)).

As we have discussed in the main text, Bell operators
$\hat{I}_{K}(\phi_{0},\theta_{0},\phi_{1},\theta_{1})$ are functions of four
angles for a given state. For $\hat{I}_{4}$ we fix the three angles
$\phi_{0},\theta_{0},\phi_{1}$ to their optimal value while leaving
$\theta=\theta_{1}$ to be variable, and show the relative quantum violation as
a function of $\theta$ in Fig 7. For the same sates and optimized measurement
directions, we plot in the same figure the violation achievable by $I_{2}$.
Clearly, we can see that the PIBI $I_{4}$ gives a relative quantum violation
of $Q_{V}^{N}/\beta_{C}^{N}=-0.1390$ for a state where $I_{2}$ fails to reveal
Bell correlations.

Figure 7: Relative quantum violation of the 2nd-order PIBI $I_{2}$ (blue) and
4th-order PIBI $I_{4}$ (red) for $N=50$ eigenstate corresponding to minimum
eigenvalue of $\hat{I}_{4}$ operator. For a given state, Bell operators
$\hat{I}_{K}(\phi_{0},\theta_{0},\phi_{1},\theta_{1})$ are associated with
four measurement angles. Here, we restricted three angles of them
$\phi_{0},\theta_{0},\phi_{1}$ to be optimal ones corresponding to minimum
relative quantum violation, and plot $Q_{V}^{N}/\beta_{C}^{N}$ changes with
the measurement angle $\theta=\theta_{1}$.