Bell inequalities tailored to the Greenberger–Horne–Zeilinger states of arbitrary local dimension

To make our considerations more accessible, we first present our construction in the case of three observers (N = 3). As already mentioned, we now denote the parties A, B and C, whereas their outcomes and measurements choices as $a,b,c$ and $x,y,z$, respectively. The departure point of our considerations is the following generalization of the Bell expression in equation (10),

$$\begin{eqnarray}&&{I}_{3,m,d}=\displaystyle \sum _{n=0}^{\lfloor d/2\rfloor -1}\left({\alpha }_{n}{{\mathbb{P}}}_{n}-{\beta }_{n}{{\mathbb{Q}}}_{n}\right).\end{eqnarray} \tag{ 25 }$$

Here, the variables defined in equations (13) and (14) simplify to

$$\begin{eqnarray}&&{X}_{\alpha ,\beta }={A}_{\alpha }-{B}_{\alpha +\beta -1}+{C}_{\beta },\,{\bar{X}}_{\alpha ,\beta }=-{A}_{\alpha +1}+{B}_{\alpha +\beta -1}-{C}_{\beta }\end{eqnarray} \tag{ 26 }$$

and therefore ${{\mathbb{P}}}_{n}$ and ${{\mathbb{Q}}}_{n}$ can be written as

$$\begin{eqnarray}&&{{\mathbb{P}}}_{n}=\displaystyle \sum _{\alpha ,\beta =1}^{m}\left[P({A}_{\alpha }-{B}_{\alpha +\beta -1}+{C}_{\beta }=n)+P({B}_{\alpha +\beta -1}-{A}_{\alpha +1}-{C}_{\beta }=n)\right]\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}&&{{\mathbb{Q}}}_{n}=\displaystyle \sum _{\alpha ,\beta =1}^{m}\left[P({A}_{\alpha }-{B}_{\alpha +\beta -1}+{C}_{\beta }=-n-1)+P({B}_{\alpha +\beta -1}-{A}_{\alpha +1}-{C}_{\beta }=-n-1)\right].\end{eqnarray} \tag{ 28 }$$

Moreover, ${\alpha }_{n}$ and ${\beta }_{n}$ are our free parameters that we are going to determine. Notice that for ${\alpha }_{n}={\beta }_{n}=[1-2n/(d-1)]$, one recovers the multipartite Bell inequalities (10). Notice also that the reason to consider the above generalization of (10), in which the free parameters multiply ${{\mathbb{P}}}_{n}$ and ${{\mathbb{Q}}}_{n}$ stems from the fact that for the state $| {\mathrm{GHZ}}_{N,d}\rangle$ and the optimal measurements (15) and (16), which we later use to find our inequalities, all probabilities contributing to ${{\mathbb{P}}}_{n}$ and ${{\mathbb{Q}}}_{n}$ are equal (see appendix for the proof).

We now want to exploit these degrees of freedom in order to construct Bell inequalities maximally violated by $| {\mathrm{GHZ}}_{3,d}\rangle$. However, in order to fully appreciate the symmetries inherent in Bell inequalities and thus significantly simplify our considerations we will write the Bell expression (25) in terms of complex correlators (20) instead of probabilities. After some short algebra one finds that

$$\begin{eqnarray}&&P({A}_{\alpha }-{B}_{\alpha +\beta -1}+{C}_{\beta }=\xi )=\displaystyle \frac{1}{d}\displaystyle \sum _{k=0}^{d-1}{\omega }^{-k\xi }\left\langle {{\mathscr{A}}}_{\alpha }^{k}\,{{\mathscr{B}}}_{\alpha +\beta -1}^{-k}\,{{\mathscr{C}}}_{\beta }^{k}\right\rangle \end{eqnarray} \tag{ 29 }$$

and

$$\begin{eqnarray}&&P(-{A}_{\alpha +1}+{B}_{\alpha +\beta -1}-{C}_{\beta }=\xi )=\displaystyle \frac{1}{d}\displaystyle \sum _{k=0}^{d-1}{\omega }^{k\xi }\left\langle {{\mathscr{A}}}_{\alpha +1}^{k}\,{{\mathscr{B}}}_{\alpha +\beta -1}^{-k}\,{{\mathscr{C}}}_{\beta }^{k}\right\rangle .\end{eqnarray} \tag{ 30 }$$

With the aid of these formulas our Bell expression (25) can be conveniently rewritten as

$$\begin{eqnarray}&&{I}_{3,m,d}=\displaystyle \frac{1}{d}\displaystyle \sum _{\alpha ,\beta =1}^{m}\displaystyle \sum _{k=0}^{d-1}\left\langle {\overline{{\mathscr{A}}}}_{\alpha }^{(k)}\,{{\mathscr{B}}}_{\alpha +\beta -1}^{-k}\,{{\mathscr{C}}}_{\beta }^{k}\right\rangle ,\end{eqnarray} \tag{ 31 }$$

where the new variables ${\overline{{\mathscr{A}}}}_{\alpha }^{(k)}$ are defined as

$$\begin{eqnarray}&&{\overline{{\mathscr{A}}}}_{\alpha }^{(k)}={a}_{k}\,{{\mathscr{A}}}_{\alpha }^{k}+{a}_{k}^{* }{{\mathscr{A}}}_{\alpha +1}^{k}\end{eqnarray} \tag{ 32 }$$

with

$$\begin{eqnarray}&&{a}_{k}=\displaystyle \sum _{n=0}^{\lfloor d/2\rfloor -1}\left[{\alpha }_{n}{\omega }^{-{kn}}-{\beta }_{n}{\omega }^{k(n+1)}\right].\end{eqnarray} \tag{ 33 }$$

Notice here that as, due to the convention, ${{\mathscr{A}}}_{m+1}=\omega {{\mathscr{A}}}_{1}$, and therefore in the particular case of $\alpha =m$, equation (32) reads ${\overline{{\mathscr{A}}}}_{m}^{(k)}={a}_{k}\,{{\mathscr{A}}}_{m}^{k}+{a}_{k}^{* }\,{\omega }^{k}{{\mathscr{A}}}_{1}^{k}$. Let us also notice that the term in equation (31) corresponding to $k=0$ is a constant and therefore it is not included in the Bell expression.

Now, to fix our free parameters ${\alpha }_{n}$ and ${\beta }_{n}$ $(n=0,\,\ldots ,\,\lfloor d/2\rfloor -1)$ we require that for the optimal observables (16) the following conditions

$$\begin{eqnarray}&&{\overline{{\mathscr{A}}}}_{\alpha }^{(k)}\otimes {{\mathscr{B}}}_{\alpha +\beta -1}^{-k}\otimes {{\mathscr{C}}}_{\beta }^{k}| {\mathrm{GHZ}}_{3,d}\rangle =| {\mathrm{GHZ}}_{3,d}\rangle \end{eqnarray} \tag{ 34 }$$

are satisfied for all $\alpha ,\beta =1,\,\ldots ,\,m$ and $k=1,\,\ldots ,\,d-1$. In other words, we want to find such ${\alpha }_{n}$ and ${\beta }_{n}$ that the resulting operator ${\overline{{\mathscr{A}}}}_{\alpha }^{(k)}\otimes {{\mathscr{B}}}_{\alpha +\beta -1}^{-k}\otimes {{\mathscr{C}}}_{\beta }^{k}$ stabilizes the GHZ state $| {\mathrm{GHZ}}_{3,d}\rangle$, or that the GHZ state is its eigenstate with eigenvalue one. To solve the above equations we need the explicit forms of the kth powers of the measurements (16). After simple algebra one finds that

$$\begin{eqnarray}&&{{\mathscr{A}}}_{x}^{k}={\omega }^{-(d-k){\gamma }_{m}(x)}\displaystyle \sum _{n=0}^{k-1}| d-k+n\rangle \langle n| +{\omega }^{k{\gamma }_{m}(x)}\displaystyle \sum _{n=k}^{d-1}| n-k\rangle \langle n| ,\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{{\mathscr{B}}}_{y}^{-k}={\left({{\mathscr{B}}}_{y}^{k}\right)}^{\dagger }={\omega }^{(d-k){\zeta }_{m}(y)}\displaystyle \sum _{n=0}^{k-1}| d-k+n\rangle \langle n| +{\omega }^{k{\zeta }_{m}(y)}\displaystyle \sum _{n=k}^{d-1}| n-k\rangle \langle n| ,\end{eqnarray} \tag{ 36 }$$

and

$$\begin{eqnarray}&&{{\mathscr{C}}}_{z}^{k}={\omega }^{-(d-k){\theta }_{m}(z)}\displaystyle \sum _{n=0}^{k-1}| d-k+n\rangle \langle n| +{\omega }^{k{\theta }_{m}(z)}\displaystyle \sum _{n=k}^{d-1}| n-k\rangle \langle n| .\end{eqnarray} \tag{ 37 }$$

By plugging (35)–(37) into the conditions (34), one obtains the following system of linear equations for the coefficients a_k:

$$\begin{eqnarray}&&\left\{\begin{array}{l}{a}_{k}{\omega }^{-k/2m}+{a}_{k}^{* }{\omega }^{k/2m}=1,\\ {a}_{k}{\omega }^{(d-k)/2m}+{a}_{k}^{* }{\omega }^{-(d-k)/2m}=1,\end{array}\right.\end{eqnarray} \tag{ 38 }$$

with $k=1,\,\ldots ,\,\lfloor d/2\rfloor$. This system can be directly solved, giving

$$\begin{eqnarray}&&{a}_{k}=\displaystyle \frac{{\omega }^{\tfrac{2k-d}{4m}}}{2\cos (\pi /2m)}\ (k=1,\,\ldots ,\,\lfloor d/2\rfloor ).\end{eqnarray} \tag{ 39 }$$

Having explicit form of a_k, we can now excavate the coefficients ${\alpha }_{n}$ and ${\beta }_{n}$ from the system of equations (33). The latter consists of $\lfloor d/2\rfloor$ equations involving $2\lfloor d/2\rfloor$ variables, meaning that it does not have a unique solution, and the solution will not be real in general. To avoid this problem we equip (33) with $\lfloor d/2\rfloor$ additional equations of the form

$$\begin{eqnarray}&&\displaystyle \sum _{n=0}^{\lfloor d/2\rfloor -1}\left({\alpha }_{n}{\omega }^{{kn}}-{\beta }_{n}{\omega }^{-(n+1)k}\right)={a}_{k}^{* }\end{eqnarray} \tag{ 40 }$$

which will enforce ${\alpha }_{n}$ and ${\beta }_{n}$ to be real. Both sets of equations (33) and (41) can be wrapped up into a single set of the form

$$\begin{eqnarray}&&\displaystyle \sum _{n=0}^{\lfloor d/2\rfloor -1}\left({\alpha }_{n}{\omega }^{-{kn}}-{\beta }_{n}{\omega }^{k(n+1)}\right)={c}_{k},\end{eqnarray} \tag{ 41 }$$

in which ${c}_{k}={a}_{k}$ for $k=1,\,\ldots ,\,\lfloor d/2\rfloor$ and ${c}_{k}={c}_{-k}^{* }$ for $k=-\lfloor d/2\rfloor ,\,\ldots ,\,-1$. In order to solve this system we consider the cases of odd and even d separately. For odd d, (41) has the same number of equations and variables, and thus we expect a unique solution, which after rather straightforward but tedious calculations is found to be

$$\begin{eqnarray}&&{\alpha }_{n}=\displaystyle \frac{1}{2d}\tan \left(\displaystyle \frac{\pi }{2m}\right)\left[{g}_{m}(n)-g(\lfloor d/2\rfloor )\right]\end{eqnarray} \tag{ 42 }$$

and

$$\begin{eqnarray}&&{\beta }_{n}=\displaystyle \frac{1}{2d}\tan \left(\displaystyle \frac{\pi }{2m}\right)\left[{g}_{m}(n+1-1/m)+{g}_{m}(\lfloor d/2\rfloor )\right]\end{eqnarray} \tag{ 43 }$$

with $n=1,\,\ldots ,\,\lfloor d/2\rfloor$ and ${g}_{m}(x):= \cot \{\pi [x+1/(2m)]/d\}$.

Then, in the case of even d the equations for $k=d/2$ and $k=-d/2$ are the same, and therefore the system (41) consists of $d-1$ equations for d variables. The additional variable we then fix in such a way that the obtained ${\alpha }_{n}$ and ${\beta }_{n}$ assume the same form as for odd d.

It turns out that the above considerations remain valid if one considers an arbitrary number of parties N.

Let us consider the same Bell expression as in (25), i.e.

$$\begin{eqnarray}&&{I}_{N,m,d}=\displaystyle \sum _{n=0}^{\lfloor d/2\rfloor -1}({\alpha }_{n}{{\mathbb{P}}}_{n}-{\beta }_{n}{{\mathbb{Q}}}_{n}),\end{eqnarray} \tag{ 44 }$$

in which now ${{\mathbb{P}}}_{n}$ and ${{\mathbb{Q}}}_{n}$ are defined for any N in (11) and (12). By using the discrete Fourier transform we can rewrite it in terms of the complex correlators as

$$\begin{eqnarray}&&{\widetilde{I}}_{N,m,d}=\displaystyle \frac{1}{d}\displaystyle \sum _{{\alpha }_{1},\,\ldots ,\,{\alpha }_{N-1}=1}^{m}\displaystyle \sum _{k=1}^{d-1}\left\langle {\overline{{\mathscr{A}}}}_{{\alpha }_{1}}^{(k)}\displaystyle \prod _{i=2}^{N}{\left({{\mathscr{A}}}_{i,{\alpha }_{i-1}+{\alpha }_{i}-1}\right)}^{{\left(-1\right)}^{i-1}k}\right\rangle ,\end{eqnarray} \tag{ 45 }$$

where ${\alpha }_{N}=1$. The variables ${\overline{{\mathscr{A}}}}_{{\alpha }_{1}}^{(k)}$ are, as before, combinations of ${{\mathscr{A}}}_{1,{\alpha }_{1}}^{k}$ and ${{\mathscr{A}}}_{1,{\alpha }_{1}+1}^{k}$ given by

$$\begin{eqnarray}&&{\overline{{\mathscr{A}}}}_{{\alpha }_{1}}^{(k)}={a}_{k}\,{{\mathscr{A}}}_{1,{\alpha }_{1}}^{k}+{a}_{k}^{* }\,{{\mathscr{A}}}_{1,{\alpha }_{1}+1}^{k}\end{eqnarray} \tag{ 46 }$$

for ${\alpha }_{1}=1,\,\ldots ,\,m$, where, as before, ${{\mathscr{A}}}_{1,m+1}=\omega {{\mathscr{A}}}_{\mathrm{1,1}}$. The coefficients a_k are defined in equation (33). Notice also that the Bell expression (66) does not contain the term corresponding to $k=0$ as the latter is only a constant that can always be moved to the classical bound; for this reason we changed the notation from ${I}_{N,m,d}$ to ${\widetilde{I}}_{N,m,d}$. The values ${I}_{N,m,d}$ and ${\widetilde{I}}_{N,m,d}$ are related by the following formula

$$\begin{eqnarray}&&{I}_{N,m,d}={\widetilde{I}}_{N,m,d}+\displaystyle \frac{{m}^{N-1}}{d}\displaystyle \sum _{n=0}^{\lfloor d/2\rfloor -1}\left({\alpha }_{n}-{\beta }_{n}\right),\end{eqnarray} \tag{ 47 }$$

and so for a particular choice of the coefficients ${\alpha }_{n}$ and ${\beta }_{n}$, given the value of our Bell expression in one representation, we can easily compute it in the other representation.

Now, the above form of the Bell expression suggests the conditions one needs to impose on the variables ${\alpha }_{n}$ and ${\beta }_{n}$ in order to obtain a Bell inequality maximally violated by the N-partite GHZ state $| {\mathrm{GHZ}}_{N,d}\rangle$. Namely, the following system of equations

$$\begin{eqnarray}&&{\overline{{\mathscr{A}}}}_{{\alpha }_{1}}^{(k)}\otimes \underset{i=2}{\overset{N}{\displaystyle \bigotimes }}{\left({{\mathscr{A}}}_{i,{\alpha }_{i-1}+{\alpha }_{i}-1}\right)}^{{\left(-1\right)}^{i-1}k}| {\mathrm{GHZ}}_{N,d}\rangle =| {\mathrm{GHZ}}_{N,d}\rangle ,\end{eqnarray} \tag{ 48 }$$

with the same conventions as above, should hold for any sequence of ${\alpha }_{i}$'s and k with the measurements being given in (16). After some tedious calculations one finds that this leads to the same system of equations for a_k as we obtained in the tripartite case (38); its solution is given in (39). Thus, our Bell inequalities for any number of parties are determined through the same coefficients ${\alpha }_{n}$ and ${\beta }_{n}$ as in the case $N=3$.

To summarize, in the probability representation our class of Bell inequalities is given by (44) with ${\alpha }_{n}$ and ${\beta }_{n}$ stated explicitly in equations (42) and (43), respectively, while in the correlator representation it is given by (66) with a_k being of the form (39). Moreover, for this choice of ${\alpha }_{n}$ and ${\beta }_{n}$, we have

$$\begin{eqnarray}&&S:= \displaystyle \sum _{n=0}^{\lfloor d/2\rfloor -1}\left({\alpha }_{n}-{\beta }_{n}\right)=\displaystyle \frac{1}{2}\left\{1-\tan \left(\displaystyle \frac{\pi }{2m}\right)\cot \left[\displaystyle \frac{\pi }{d}\left(\left\lfloor \displaystyle \frac{d}{2}\right\rfloor +\displaystyle \frac{1}{2m}\right)\right]\right\},\end{eqnarray} \tag{ 49 }$$

and consequently the Bell expression in the probability and correlator forms are related through the following simple formula

$$\begin{eqnarray}&&{I}_{N,m,d}={\widetilde{I}}_{N,m,d}+(2{m}^{N-1}/d)S.\end{eqnarray} \tag{ 50 }$$

To conclude this part, let us mention the Bell expressions obtained here belong to a more general family of Bell expressions considered in [26], some of which independently discovered in [27]. However, in these works it is only shown that high-dimensional GHZ states can exhibit fully random and genuinely multipartite quantum correlations. Our method is such that the inequalities are built from the property that the multiqudit GHZ state and given measurements maximally violate it.

Let us begin by noting that our Bell expression ${I}_{N,m,d}$ can be written in a simpler form as

$$\begin{eqnarray}&&{I}_{N,m,d}=\displaystyle \sum _{n=0}^{d-1}{\hat{\alpha }}_{n}{{\mathbb{P}}}_{n},\end{eqnarray} \tag{ 51 }$$

where ${{\mathbb{P}}}_{n}$ is given in equation (11) and ${\hat{\alpha }}_{n}={\alpha }_{n}$ for $n=0,\,\ldots ,\,\lfloor d/2\rfloor -1$ and ${\hat{\alpha }}_{n}=-{\beta }_{d-1-n}$ for $n=\lfloor d/2\rfloor ,\,\ldots ,\,d-1$ (notice that in the odd d case ${\alpha }_{\lfloor d/2\rfloor }={\beta }_{\lfloor d/2\rfloor }=0$).

Let us also recall that to compute the maximal classical value ${\beta }_{{ \mathcal L }}$ of our Bell expressions ${I}_{N,m,d}$ it is enough to maximize the latter over the vertices of ${{ \mathcal L }}_{N,m,d}$, or, in other words, over all deterministic assignments ${A}_{i,{x}_{i}}\in [d]$ for ${x}_{i}=1,\,\ldots ,\,m$ and $i=1,\,\ldots ,\,N$, i.e.

$$\begin{eqnarray}&&{\beta }_{{ \mathcal L }}=\mathop{\max }\limits_{{\{{A}_{i,k}\in [d]\}}_{i\in [N],k\in [M]}}{I}_{N,m,d},\end{eqnarray} \tag{ 52 }$$

for which ${I}_{N,m,d}$ given in (51) rewrites as

$$\begin{eqnarray}&&{I}_{N,m,d}=\displaystyle \sum _{n=0}^{d-1}{\hat{\alpha }}_{n}\displaystyle \sum _{{\alpha }_{1},\ldots ,{\alpha }_{N-1}=1}^{m}\left[\delta ({X}_{{\alpha }_{1},\ldots ,{\alpha }_{N-1}},n)+\delta ({\overline{X}}_{{\alpha }_{1},\ldots ,{\alpha }_{N-1}},n)\right],\end{eqnarray} \tag{ 53 }$$

where $\delta (\cdot ,\,\cdot )$ denotes the Kronecker delta, whereas the variables X and $\overline{X}$ are defined in equations (13) and (14).

To facilitate the computation of the classical value we want to express the maximum in (52) in terms of the X and $\overline{X}$ variables, instead of ${A}_{i,{x}_{i}}$. To do this, we need to remove all the linear dependencies between ${X}_{{\alpha }_{1},\ldots ,{\alpha }_{N-1}}$ and ${\overline{X}}_{{\alpha }_{1},\ldots ,{\alpha }_{N-1}}$. Let us illustrate what we mean by this with the bipartite case in which the classical value was computed analytically for any m and d in [19].

Let us then assume that $N=2$ and notice that the variables ${X}_{\alpha }$ and ${\overline{X}}_{\alpha }$ are related to ${A}_{\alpha }$ and ${B}_{\alpha }$ by the following formula

$$\begin{eqnarray}&&\left(\begin{array}{c}1\\ {X}_{1}\\ \vdots \\ {X}_{m}\\ {\overline{X}}_{1}\\ \vdots \\ {\overline{X}}_{m}\end{array}\right)=H\left(\begin{array}{c}1\\ {A}_{1}\\ \vdots \\ {A}_{m}\\ {B}_{1}\\ \vdots {B}_{m}\end{array}\right),\end{eqnarray} \tag{ 54 }$$

where

$$\begin{eqnarray}H=\left(\begin{array}{ccc}1 & 0 & 0\\ 0 & {\mathbb{1}} & -{\mathbb{1}}\\ b & -{ \mathcal X } & {\mathbb{1}}\end{array}\right),\end{eqnarray} \tag{ 55 }$$

$b={\left(0,0,\ldots ,0,-1\right)}^{T}$ and ${ \mathcal X }={\sum }_{i=0}^{m-1}| i\rangle \langle i+1|$. In order to find a linearly independent set of ${X}_{\alpha }$ and ${\overline{X}}_{\alpha }$, we want to find all linear combinations such that

$$\begin{eqnarray}&&a+\displaystyle \sum _{\alpha }\left[{g}_{\alpha }{X}_{\alpha }+{h}_{\alpha }{\overline{X}}_{\alpha }\right]=0,\end{eqnarray} \tag{ 56 }$$

where a, ${g}_{\alpha }$ and ${h}_{\alpha }$ are some coefficients to be determined. For this purpose, we want to determine the kernel of H^T, which in this case consists of a single vector ${\left(1,1,\ldots ,1\right)}^{T}$. Consequently, we arrive at the condition for the bipartite Bell expression ${I}_{2,m,d}$ constructed in [19], which is

$$\begin{eqnarray}&&\displaystyle \sum _{\alpha =0}^{m-1}\left[{X}_{\alpha }+{\overline{X}}_{\alpha }\right]\equiv -1\mathrm{mod}d.\end{eqnarray} \tag{ 57 }$$

As proven in [19], this allows to nest the optimization of the classical bound and use a dynamic programming procedure to find it efficiently. Indeed, one can re-express the optimization in (52) in terms of the ${X}_{\alpha },{\overline{X}}_{\alpha }$, but with the condition stemming from (57).

This allows us to write

$$\begin{eqnarray}&&{\beta }_{{ \mathcal L }}=\mathop{\max }\limits_{\left\{{X}_{\alpha },{\overline{X}}_{\alpha }\in [d]:\displaystyle \sum _{\alpha }{X}_{\alpha }+{\overline{X}}_{\alpha }\equiv -1\mathrm{mod}d\right\}}{I}_{N,m,d},\end{eqnarray} \tag{ 58 }$$

and eliminate the variables in the optimization successively thanks to the form of (53). The corresponding classical bound of ${I}_{2,m,d}$ is then found to be [19]:

$$\begin{eqnarray}&&{\beta }_{{ \mathcal L }}^{2,m,d}=\displaystyle \frac{1}{2d}\tan \left(\displaystyle \frac{\pi }{2m}\right)\left[(2m-1){g}_{m}(0)-{g}_{m}\left(1-\displaystyle \frac{1}{m}\right)-2{{mg}}_{m}\left(\left\lfloor \displaystyle \frac{d}{2}\right\rfloor \right)\right].\end{eqnarray} \tag{ 59 }$$

In what follows we attempt to formalize the procedure to find the classical bound of the inequality for a larger system size. Here we outline the steps, but several obstacles arise in the multipartite case that we currently do not see how to overcome.

Let us consider an arbitrary number of observers N. Likewise, we wish to find all $a,{g}_{{\boldsymbol{\alpha }}},{h}_{{\boldsymbol{\alpha }}}\in {{\mathbb{Z}}}_{d}$ such that

$$\begin{eqnarray}&&a+\displaystyle \sum _{{\boldsymbol{\alpha }}}\left[{g}_{{\boldsymbol{\alpha }}}{X}_{{\boldsymbol{\alpha }}}+{h}_{{\boldsymbol{\alpha }}}{\overline{X}}_{{\boldsymbol{\alpha }}}\right]\equiv 0\mathrm{mod}d.\end{eqnarray} \tag{ 60 }$$

Note however, that now ${\boldsymbol{\alpha }}$ is a vector and therefore one expects the solution set to be multidimensional. In other words, H is a square matrix only for $N=2$. If $N\gt 2$ we shall have on the left hand side of (54) a $(1+2{m}^{N-1})$-component vector of ${X}_{{\boldsymbol{\alpha }}}$'s and ${\overline{X}}_{{\boldsymbol{\alpha }}}$'s, whereas on the left hand side we shall have a ($1+{mN}$)-component vector of observables. Therefore, the number of terms in the kernel of H^T will grow exponentially with N, implying there will be an exponential number of non-trivial constraints in the maximization analogous to (58). In order to find them, using equations (13) and (14), the above equation can be expanded as

$$\begin{eqnarray}&&a+\displaystyle \sum _{{\boldsymbol{\alpha }}}\left[{g}_{{\boldsymbol{\alpha }}}{A}_{1,{\alpha }_{1}}-{h}_{{\boldsymbol{\alpha }}}{A}_{1,{\alpha }_{1}+1}\right]+\displaystyle \sum _{i=2}^{N}{\left(-1\right)}^{i-1}\displaystyle \sum _{{\boldsymbol{\alpha }}}\left[{g}_{{\boldsymbol{\alpha }}}-{h}_{{\boldsymbol{\alpha }}}\right]{A}_{i,{\alpha }_{i-1}+{\alpha }_{i}-1}\equiv 0\mathrm{mod}d.\end{eqnarray} \tag{ 61 }$$

This gives a set of equations that make the coefficients in front of ${A}_{N,{\alpha }_{N-1}}$ congruent to $0\,\mathrm{mod}\,d$:

$$\begin{eqnarray}&&\displaystyle \sum _{{\boldsymbol{\alpha }}:{\alpha }_{N-1}=k}\left[{g}_{{\boldsymbol{\alpha }}}-{h}_{{\boldsymbol{\alpha }}}\right]\equiv 0\mathrm{mod}d\ (1\leqslant k\leqslant m).\end{eqnarray} \tag{ 62 }$$

Similarly, for $1\lt j\lt N$, we make the coefficient in front of ${A}_{j,k}$ congruent to $0\,\mathrm{mod}\,d$:

$$\begin{eqnarray}&&\displaystyle \sum _{{\boldsymbol{\alpha }}:{\alpha }_{j-1}+{\alpha }_{j}-1=k}\left[{g}_{{\boldsymbol{\alpha }}}-{h}_{{\boldsymbol{\alpha }}}\right]\equiv 0\mathrm{mod}d\ (1\leqslant k\leqslant m).\end{eqnarray} \tag{ 63 }$$

Then, the coefficient that multiplies ${A}_{1,{\alpha }_{1}}$ is

$$\begin{eqnarray}&&\displaystyle \sum _{{\boldsymbol{\alpha }}:{\alpha }_{1}=k}\left[{g}_{{\alpha }_{1}=k,{\boldsymbol{\alpha }}^{\prime} }-{h}_{{\alpha }_{1}=k-1,{\boldsymbol{\alpha }}^{\prime} }\right]\equiv 0\mathrm{mod}d\ (1\leqslant k\leqslant m),\end{eqnarray} \tag{ 64 }$$

where ${\boldsymbol{\alpha }}^{\prime} \equiv {\alpha }_{2},\,\ldots ,\,{\alpha }_{N-1}$. We have also an equation for the constant term in equation (60). Here we have to take into consideration that ${A}_{i,m-1+k}={A}_{i,k}+1$ for $k\gt 0$ and any $i=1,\,\ldots ,\,N$.

In order to appropriately do the substitutions of the ${X}_{{\boldsymbol{\alpha }}}$'s and ${\overline{X}}_{{\boldsymbol{\alpha }}}$'s, one needs to perform Gauss elimination on a basis of this kernel. However, we note that the equations (60) are over ${{\mathbb{Z}}}_{d}$, which is a field only of d is a prime number. Therefore, for some values of d, inverses may not exist, and it can be much more complicated to obtain a good basis of $\ker {H}^{T}$. This was not a problem for $N=2$ as $\ker {H}^{T}$ was generated simply by a vector of ones. In addition, it is unclear how to later exploit the properties of the g function that were used in [19] in order to find an analytical form of $\beta {{ \mathcal L }}^{N,m,d}$ for $N\gt 2$.

Thus, even though the procedure above helps to slightly reduce the complexity of the optimization in some particular cases with a few number of particles, we resort to numerics in the general case. As pointed out in theorem 1, we have that ${\beta }_{{ \mathcal L }}^{3,2,d}={\beta }_{{ \mathcal S }}^{3,2,d}$ and its value is thus given by expression (59). For completeness, we include the classical bound values in the simplest Bell scenarios in tables 1 and 2.

Table 1.  Maximal classical and Svetlichny values of ${I}_{N,m,d}$ for $N=3$, $m=2,3$ and $d=2,3,4$.

$N=3$, $m=2$
d	2	3	4
${\beta }_{{ \mathcal L }}^{3,2,d}={\beta }_{{ \mathcal S }}^{3,2,d}$	4.2426	3.0416	3.5953
$N=3$, $m=3$
${\beta }_{{ \mathcal L }}^{3,3,d}$	$\tfrac{13}{\sqrt{3}}$	$\tfrac{1}{6\sqrt{3}}\left[13\cot \left(\tfrac{\pi }{18}\right)-17\tan \left(\tfrac{\pi }{9}\right)-4\tan \left(\tfrac{2\pi }{9}\right)\right]$	$\tfrac{-10+17\sqrt{2}+14\sqrt{6}}{4\sqrt{3}}$
	= 7.5056	= 6.1760	= 6.9765
${\beta }_{{ \mathcal S }}^{3,3,d}$	8.6603	7.3132	8.1115

Table 2.  Maximal classical and Svetlichny values of ${I}_{N,m,d}$ for $N=4$, $m=2,3$ and $d=2,3,4$. It was not possible for us to compute the value of ${\beta }_{{ \mathcal L }}^{4,3,4}$ due to its computational complexity.

$N=4$, $m=2$
d	$d=2$	$d=3$	$d=4$
${\beta }_{{ \mathcal L }}^{4,2,d}$	$\tfrac{5}{2}\left[\cot \left(\tfrac{\pi }{8}\right)+\tan \left(\tfrac{\pi }{8}\right)\right]$	$\tfrac{10}{\sqrt{3}}+\tfrac{5}{6}(-3+\sqrt{3})$	$\tfrac{1}{8}[10\cot \left(\tfrac{\pi }{16}\right)-5\cot \left(\tfrac{3\pi }{16}\right)$
			$+16\tan \left(\tfrac{\pi }{16}\right)+\tan \left(\tfrac{3\pi }{16}\right)]$
	=7.0711	=4.7169	=5.8301
${\beta }_{{ \mathcal S }}^{4,2,d}$	8.4853	6.0829	7.1905
$N=4$, $m=3$
${\beta }_{{ \mathcal L }}^{4,3,d}$	$\tfrac{35}{\sqrt{3}}$	$\tfrac{7}{6\sqrt{3}}\left[5\cot \left(\tfrac{\pi }{18}\right)-7\tan \left(\tfrac{\pi }{9}\right)-2\tan \left(\tfrac{2\pi }{9}\right)\right]$	(n.a.)
	=20.2073	=16.2537
${\beta }_{{ \mathcal S }}^{4,3,d}$	25.9808	21.9394	24.3345

On the other hand, it is quite direct to obtain an upper bound on the maximal value of ${I}_{N,m,d}$ over the Svetlichny correlations. Precisely, the results of [26] allow us to state the following theorem.

Theorem 1. The Svetlichny bound of ${I}_{N,m,d}$ is bounded from above as ${\beta }_{{ \mathcal S }}^{N,m,d}\leqslant {m}^{N-2}{\beta }_{{ \mathcal L }}^{2,m,d}$, where ${\beta }_{{ \mathcal L }}^{2,m,d}$ is the classical bound of the bipartite Bell inequality given explicitly in (59).

Proof. The proof is given in [26]. □

It is worth mentioning that for the case $N=3$ and $m=2$ the bound ${\beta }_{{ \mathcal S }}^{3,2,d}$ is also saturated by fully product probability distribution (see also table 1)

$$\begin{eqnarray}&&{p}_{\mathrm{loc}}({\bf{a}}| {\bf{x}})={p}_{A}(a| x){p}_{B}(b| y){p}_{C}(c| z)\end{eqnarray} \tag{ 65 }$$

such that ${p}_{A}(0| x)={p}_{A}(0| y)={p}_{C}(0| z)=1$ for all $x,y,z$. So, ${\beta }_{{ \mathcal S }}^{3,2,d}$ is also the classical bound of the corresponding Bell inequality. In general, however, the classical and Svetlichny bounds differ.

Let us now move on to the quantum and nonsignaling bounds.

Theorem 2. The maximal quantum value of ${\widetilde{I}}_{N,m,d}$ is ${\widetilde{\beta }}_{Q}^{N,m,d}={m}^{N-1}(d-1)/d$.

Proof. To prove that ${\widetilde{\beta }}_{Q}^{N,m,d}$ is an upper bound on the maximal quantum value of ${I}_{N,m,d}$ we can follow the method of [26]. Here, however, we follow an alternative approach exploiting the sum-of-squares decompositon of the shifted Bell operator, which might be of use for such applications of nonlocality as self-testing.

To this end, let us consider a Bell operator ${{ \mathcal B }}_{N,m,d}$ constructed from the Bell expression ${\widetilde{I}}_{N,d,m}$ with some observables ${{\mathscr{A}}}_{i,{x}_{i}}$ (that is, unitary operators such that ${{\mathscr{A}}}_{i,{x}_{i}}^{d}={\mathbb{1}}$):

$$\begin{eqnarray}&&{{B}}_{N,m,d}\,:=\displaystyle \frac{1}{d}\displaystyle \sum _{{\alpha }_{1},\,\ldots ,\,{\alpha }_{N-1}=1}^{m}\displaystyle \sum _{k=1}^{d-1}{\overline{{\mathscr{A}}}}_{{\alpha }_{1}}^{(k)}\otimes \mathop{\mathop{\displaystyle \bigotimes }\limits^{N}}\limits_{i=2}{\left({{\mathscr{A}}}_{i,{\alpha }_{i-1}+{\alpha }_{i}-1}\right)}^{{\left(-1\right)}^{i-1}k}.\end{eqnarray} \tag{ 66 }$$

Our aim now is to prove that the operator ${\widetilde{\beta }}_{{ \mathcal Q }}^{N,m,d}{\mathbb{1}}-{{ \mathcal B }}_{N,m,d}$ is positive semi-definite for arbitrary observables ${{\mathscr{A}}}_{i,{x}_{i}}$ with the identity operator ${\mathbb{1}}$ being of defined on the corresponding Hilbert space. This can be achieved by decomposing ${\widetilde{\beta }}_{{ \mathcal Q }}^{N,m,d}{\mathbb{1}}-{{ \mathcal B }}_{N,m,d}$ into a sum of squares. To be more precise, let us first consider the simpler case of $m=2$ and introduce the following operators

$$\begin{eqnarray}&&{P}_{{\alpha }_{1},\,\ldots ,\,{\alpha }_{N-1}}^{(k)}={\mathbb{1}}-{\overline{{\mathscr{A}}}}_{{\alpha }_{1}}^{(k)}\otimes \underset{i=2}{\overset{N}{\displaystyle \bigotimes }}{\left({{\mathscr{A}}}_{i,{\alpha }_{i-1}+{\alpha }_{i}-1}\right)}^{{\left(-1\right)}^{i-1}k}.\end{eqnarray} \tag{ 67 }$$

Then, by a direct check one finds that the following decomposition

$$\begin{eqnarray}&&{\widetilde{\beta }}_{{ \mathcal Q }}^{N,2,d}{\mathbb{1}}-{{ \mathcal B }}_{N,2,d}=\displaystyle \frac{1}{2d}\displaystyle \sum _{{\alpha }_{1},\,\ldots ,\,{\alpha }_{N-1}=1}^{m}\displaystyle \sum _{k=1}^{d-1}{\left[{P}_{{\alpha }_{1},\ldots ,{\alpha }_{N-1}}^{(k)}\right]}^{\dagger }{P}_{{\alpha }_{1},\,\ldots ,\,{\alpha }_{N-1}}^{(k)}\end{eqnarray} \tag{ 68 }$$

holds true. In the case of arbitrary number of measurements, the above sum of squares needs to be slightly modified. Let us introduce the following operators

$$\begin{eqnarray}&&{T}_{\alpha }^{(k)}={\mu }_{\alpha ,k}^{* }\,{{\mathscr{A}}}_{1,2}^{k}+{\nu }_{\alpha ,k}^{* }\,{{\mathscr{A}}}_{1,\alpha +2}^{k}+{\tau }_{\alpha ,k}\,{{\mathscr{A}}}_{1,\alpha +3}^{k}\end{eqnarray} \tag{ 69 }$$

for $\alpha =1,\,\ldots ,\,m-2$ and $k=1,\,\ldots ,\,d-1$, where the coefficients are defined as

$$\begin{eqnarray}\begin{array}{rcl}{\mu }_{\alpha ,k} & = & \displaystyle \frac{{\omega }^{(\alpha +1)(d-2k)/2m}}{2\cos (\pi /2m)}\displaystyle \frac{\sin (\pi /m)}{\sqrt{\sin (\pi \alpha /m)\sin \left[\pi (\alpha +1)/m\right]}},\\ {\nu }_{\alpha ,k} & = & -\displaystyle \frac{{\omega }^{(d-2k)/2m}}{2\cos (\pi /2m)}\sqrt{\displaystyle \frac{\sin \left[\pi (\alpha +1)/m\right]}{\sin (\pi \alpha /m)}},\\ {\tau }_{\alpha ,k} & = & \displaystyle \frac{1}{2\cos (\pi /2m)}\sqrt{\displaystyle \frac{\sin (\pi \alpha /m)}{\sin \left[\pi (\alpha +1)/m\right]}}=-\displaystyle \frac{{\omega }^{(d-2k)/2m}}{4{\cos }^{2}(\pi /2m)}{\nu }_{\alpha ,k}^{-1},\end{array}\end{eqnarray} \tag{ 70 }$$

for $i=1,\ldots ,m-3$ and $k=1,\ldots ,d-1$, while for $i=m-2$ and $k=1,\,\ldots ,\,d-1$ they are given by

$$\begin{eqnarray}\begin{array}{rcl}{\mu }_{m-2,k} & = & -\displaystyle \frac{{\omega }^{-(d-2k)/2m}}{2\sqrt{2}\cos (\pi /2m)\sqrt{\cos (\pi /m)}},\\ {\nu }_{m-2,k} & = & -\displaystyle \frac{{\omega }^{k}{\omega }^{(d-2k)/2m}}{2\sqrt{2}\cos (\pi /2m)\sqrt{\cos (\pi /m)}},\\ {\tau }_{m-2,k} & = & \displaystyle \frac{\sqrt{\cos (\pi /m)}}{\sqrt{2}\cos (\pi /2m)}.\end{array}\end{eqnarray} \tag{ 71 }$$

Then, the sum of squares is given by

$$\begin{eqnarray}&&{\beta }_{{ \mathcal Q }}^{N,m,d}{\mathbb{1}}-{{ \mathcal B }}_{N,m,d}=\displaystyle \frac{1}{2d}\displaystyle \sum _{{\alpha }_{1},\,\ldots ,\,{\alpha }_{N-1}=1}^{m}\displaystyle \sum _{k=1}^{d-1}{\left[{P}_{{\alpha }_{1},\ldots ,{\alpha }_{N-1}}^{(k)}\right]}^{\dagger }{P}_{{\alpha }_{1},\,\ldots ,\,{\alpha }_{N-1}}^{(k)}+\displaystyle \frac{{m}^{N-2}}{2d}\displaystyle \sum _{\alpha =1}^{m-2}\displaystyle \sum _{k=1}^{d-1}{\left[{T}_{\alpha }^{(k)}\right]}^{\dagger }{T}_{\alpha }^{(k)}.\end{eqnarray} \tag{ 72 }$$

To conclude the proof, let us notice that for the state $| {\mathrm{GHZ}}_{N,d}\rangle$ and the measurements (15) and (16) the value of ${\widetilde{I}}_{N,m,d}$ is clearly ${m}^{N-1}(d-1)/d$, which follows from the fact that for this realisation each correlator in ${\widetilde{I}}_{N,m,d}$ assumes value one (see equation (48)). ${\beta }_{{ \mathcal Q }}^{N,m,d}$ is thus the maximal quantum value of ${\widetilde{I}}_{N,m,d}$.□

Theorem 3. The maximal nonsignaling value of ${I}_{N,m,d}$ equals its algebraic bound and it is given by ${\beta }_{{ \mathcal N }}^{N,m,d}=2{m}^{N-1}{\alpha }_{0}$.

Proof. To prove this statement we use the form of ${I}_{N,m,d}$ given in equation (51). As shown in [19] (see the supplemental material), ${\alpha }_{0}\geqslant {\alpha }_{n}$ for any $0\leqslant n\leqslant d-1$, and consequently one obtains the following bound by putting all the terms in ${{\mathbb{P}}}_{0}$ equal to one:

$$\begin{eqnarray}&&{I}_{N,m,d}\leqslant 2{m}^{N-1}{\alpha }_{0}.\end{eqnarray} \tag{ 73 }$$

Now, there exists a nonsignaling probability distribution for which this inequality is saturated. For the first portion of measurement choices it is defined as

$$\begin{eqnarray}&&\begin{array}{l}p({a}_{1},\,\ldots ,\,{a}_{N}| {\alpha }_{1},{\alpha }_{1}+{\alpha }_{2}-1,\,\ldots ,\,{\alpha }_{N-2}+{\alpha }_{N-1}-1,{\alpha }_{N-1})\\ \ \ \ =\ \left\{\begin{array}{lc}\displaystyle \frac{1}{{d}^{N-1}}, & \displaystyle \sum _{i=1}^{N}{\left(-1\right)}^{i-1}{a}_{i}=f({\alpha }_{1},\,\ldots ,\,{\alpha }_{N-1})\\ 0, & \mathrm{otherwise}\end{array}\right.,\end{array}\end{eqnarray} \tag{ 74 }$$

for ${\alpha }_{i}=0,\,\ldots ,\,m-1$ and $i=1,\,\ldots ,\,N-1$, where the function f is defined as

$$\begin{eqnarray}&&f({\alpha }_{1},\ldots ,{\alpha }_{N-1})=\displaystyle \sum _{i=1}^{N-2}{\left(-1\right)}^{i-1}H({\alpha }_{i}+{\alpha }_{i+1}-m-1).\end{eqnarray} \tag{ 75 }$$

with H being the discrete Heaviside step function, defined as $H(x)=1$ if $x\geqslant 0$ and $H(x)=0$ otherwise. This function f is introduced to take into account the convention ${A}_{i,m+k}={A}_{i,k}+1$, which modifies the condition defining the probabilities in the Bell expression. Indeed, looking at the expression (11), one sees that if for all $i=1,\,\ldots ,\,N-2$, ${\alpha }_{i}+{\alpha }_{i+1}-1\leqslant m-1$, then $f=0$, but if for some j's, ${\alpha }_{j}+{\alpha }_{j+1}-1\gt m-1$, then f could be different than 0.

Then, for the other portion of measurement choices it is defined as

$$\begin{eqnarray}&&\begin{array}{l}p({a}_{1},\,\ldots ,\,{a}_{N}| {\alpha }_{1}+1,{\alpha }_{1}+{\alpha }_{2}-1,\,\ldots ,\,{\alpha }_{N-2}+{\alpha }_{N-1}-1,{\alpha }_{N-1})\\ \ \ \ =\ \left\{\begin{array}{lc}\displaystyle \frac{1}{{d}^{N-1}}, & \displaystyle \sum _{i=1}^{N}{\left(-1\right)}^{i-1}{a}_{i}=\widetilde{f}({\alpha }_{1},\,\ldots ,\,{\alpha }_{N-1})\\ 0, & \mathrm{otherwise}\end{array}\right.,\end{array}\end{eqnarray} \tag{ 76 }$$

where the function $\widetilde{f}$ is defined in the same way as f, but also takes into account that ${\alpha }_{1}+1$ can be larger than $m-1$. Thus

$$\begin{eqnarray}&&\widetilde{f}({\alpha }_{1},\ldots ,{\alpha }_{N-1})=-H({\alpha }_{1}+1-m)+f({\alpha }_{1},\ldots ,{\alpha }_{N-1}).\end{eqnarray} \tag{ 77 }$$

For all the remaining choices of measurements we define

$$\begin{eqnarray}&&p({a}_{1},\,\ldots ,\,{a}_{N}| {\alpha }_{1},\,\ldots ,\,{\alpha }_{N})=\displaystyle \frac{1}{{d}^{N}}.\end{eqnarray} \tag{ 78 }$$

Let us now recall the no-signalling principle for many parties. For the distribution of elements $p({a}_{1},\ldots ,{a}_{N}| {x}_{1},\ldots ,{x}_{N})$, the marginal $p({a}_{{i}_{1}},\ldots ,{a}_{{i}_{k}}| {x}_{{i}_{1}},\ldots ,{x}_{{i}_{k}})$ for any subset $\{{i}_{1},\ldots ,{i}_{k}\}$ of the N parties should be independent of the measurement settings of the remaining N − k parties:

$$\begin{eqnarray}&&p({a}_{{i}_{1}},\ldots ,{a}_{{i}_{k}}| {x}_{1},\ldots ,{x}_{N})=p({a}_{{i}_{1}},\ldots ,{a}_{{i}_{k}}| {x}_{{i}_{1}},\ldots ,{x}_{{i}_{k}}).\end{eqnarray} \tag{ 79 }$$

It is not difficult to verify that the distribution presented above obeys the no-signaling principle. Tracing out a single subsystem one always obtains a maximally random probability distribution.□

Let us notice that in theorem 2 we compute the maximal quantum value of our Bell expression in the correlator representation (66), whereas in theorem 3 we compute the maximal nonsignalling value in the probability picture. To obtain these values in the other picture one can use equation (50).

Let us also notice that both these values are related to the same values of the bipartite SATWAP Bell inequality by a factor ${m}^{N-2}$. Thus, the results about the relative scaling of these bounds from [19] holds for many parties. In particular ${\mathrm{lim}}_{d\to \infty }{\beta }_{{ \mathcal N }{ \mathcal S }}^{N,m,d}/{\beta }_{{ \mathcal Q }}^{N,m,d}=(2m/\pi )\tan (\pi /2m)$, and so the separation between the maximal nonsignaling and quantum values becomes smaller for larger m.

On the other hand, the classical value ${\beta }_{{ \mathcal L }}^{N,m,d}$ does not seem to obey ${\beta }_{{ \mathcal L }}^{N,m,d}={m}^{N-2}{\beta }_{{ \mathcal L }}^{2,m,d}$ (see tables 1 and 2), and therefore one can expect that the behaviour of the ratio ${\beta }_{{ \mathcal Q }}^{N,m,d}/{\beta }_{{ \mathcal L }}^{N,m,d}$ will exhibit a behaviour for large d or m different from ${\beta }_{{ \mathcal Q }}^{2,m,d}/{\beta }_{{ \mathcal L }}^{2,m,d}$ (see [19]).

Let us here briefly discuss the form of our Bell expressions in the special cases of $d=2$ and any $m\geqslant 2$, and $m=2$ and any $d\geqslant 2$. In the first one, equation (44) simplifies to

$$\begin{eqnarray}&&{I}_{N,m,2}={\alpha }_{0}{{\mathbb{P}}}_{0},\end{eqnarray} \tag{ 80 }$$

where ${\alpha }_{0}=1/[2\cos (\pi /2m)]$ (notice also that ${\beta }_{0}=0$). Then, in the correlator picture there is a single number ${a}_{1}={\alpha }_{0}=1/[2\cos (\pi /2m)]$ and therefore

$$\begin{eqnarray}&&{\overline{{\mathscr{A}}}}_{{\alpha }_{1}}^{(1)}={a}_{1}({{\mathscr{A}}}_{{\alpha }_{1}}+{{\mathscr{A}}}_{{\alpha }_{1}+1}),\end{eqnarray} \tag{ 81 }$$

where ${{\mathscr{A}}}_{1,m+1}=-{{\mathscr{A}}}_{\mathrm{1,1}}$. Then, the Bell inequality in the correlator picture can be written as

$$\begin{eqnarray}&&{\widetilde{I}}_{N,m,2}=\displaystyle \frac{{a}_{1}}{2}\displaystyle \sum _{{\alpha }_{1},\,\ldots ,\,{\alpha }_{N-1}=1}^{m}\left[\left\langle {{\mathscr{A}}}_{1,{\alpha }_{1}}\displaystyle \prod _{i=2}^{N}{{\mathscr{A}}}_{i,{\alpha }_{i-1}+{\alpha }_{i}-1}\right\rangle +\left\langle {{\mathscr{A}}}_{1,{\alpha }_{1}+1}\displaystyle \prod _{i=2}^{N}{{\mathscr{A}}}_{i,{\alpha }_{i-1}+{\alpha }_{i}-1}\right\rangle \right]\leqslant {\beta }_{{ \mathcal L }}^{N,m,2},\end{eqnarray} \tag{ 82 }$$

where ${\alpha }_{N}=1$. This a generalization of the bipartite chained Bell inequalities [33] to the multipartite scenario (see also [27] for an extension in a similar spirit, albeit in which the GHZ state does not yield the maximal violation in general). In fact, for $N=2$, after dividing by ${a}_{1}/2$, one obtains

$$\begin{eqnarray}&&{\widetilde{I}}_{2,m,2}=\displaystyle \sum _{\alpha =1}^{m}\left[\langle {{\mathscr{A}}}_{\alpha }{{\mathscr{B}}}_{\alpha }\rangle +\langle {{\mathscr{A}}}_{\alpha +1}{{\mathscr{B}}}_{\alpha }\rangle \right]\leqslant 2(m-1),\end{eqnarray} \tag{ 83 }$$

where ${{\mathscr{A}}}_{m+1}=-{{\mathscr{A}}}_{1}$, which is the chained Bell inequality [33].

In the case of $m=2$ and any d, ${I}_{N,2,d}$ is given in equation (44) with the coefficients in the probability picture simplifying to

$$\begin{eqnarray}&&{\alpha }_{k}=\displaystyle \frac{1}{2d}\left[{g}_{2}(k)+{\left(-1\right)}^{d}\tan \left(\displaystyle \frac{\pi }{4d}\right)\right],\ {\beta }_{k}=\displaystyle \frac{1}{2d}\left[{g}_{2}(k+1/2)-{\left(-1\right)}^{d}\tan \left(\displaystyle \frac{\pi }{4d}\right)\right],\end{eqnarray} \tag{ 84 }$$

whereas in the correlator picture to ${a}_{k}={\omega }^{(2k-8)/d}/\sqrt{2}$.

In this special case, equation (44) gives a class of Bell inequalities involving two parameters ${\alpha }_{0}{{\mathbb{P}}}_{0}-{\beta }_{0}{{\mathbb{Q}}}_{0}\leqslant C$, where C is the maximal classical value. However, we can always divide the whole expression by one of them, say ${\alpha }_{0}$ (provided that it is positive), reducing the number of free parameters to one. As a result we obtain the following class of Bell inequalities

$$\begin{eqnarray}&&\begin{array}{l}{J}_{\mathrm{2,2,3}}(\xi ):= P({A}_{1}={B}_{1})+P({A}_{2}={B}_{2})+P({A}_{1}={B}_{2}-1)+P({A}_{2}={B}_{1})\\ \ \ -\ \xi [P({A}_{1}={B}_{1}-1)+P({A}_{2}={B}_{2}-1)+P({A}_{1}={B}_{2})+P({A}_{2}={B}_{1}+1)]\leqslant {C}_{3}(\xi ),\end{array}\end{eqnarray} \tag{ 85 }$$

parametrized by a single parameter ξ, defined as $\xi ={\beta }_{0}/{\alpha }_{0}$. It turns out that the classical bound of these inequalities can be easily found by maximizing ${J}_{\mathrm{2,2,3}}(\xi$) over all local deterministic strategies, which gives

$$\begin{eqnarray}{C}_{3}(\xi )=\left\{\begin{array}{cl}-4\xi , & \mathrm{if}\,\,\xi \leqslant -1,\\ 3-\xi , & \mathrm{if}\,\,-1\leqslant \xi \leqslant 1,\\ 2, & \mathrm{if}\,\,\xi \geqslant 1.\end{array}\right.\end{eqnarray} \tag{ 86 }$$

Moreover, numerical tests using the Navascués–Pironio–Acín (NPA) hierarchy [34] indicate that for $\xi \leqslant -1$, the Bell inequality (85) is trivial, meaning that its maximal quantum violation equals its classical bound. Consequently, in what follows we will concentrate on the case $\xi \gt -1$. It is not difficult to see that for $\xi =1$ the class (85) reproduces the well-known CGLMP Bell inequality [8], which is known to be maximally violated by the partially entangled state [35]:

$$\begin{eqnarray}&&| {\psi }_{\gamma }\rangle =\displaystyle \frac{1}{\sqrt{2+{\gamma }^{2}}}(| 00\rangle +\gamma | 11\rangle +| 22\rangle )\end{eqnarray} \tag{ 87 }$$

with $\gamma =(\sqrt{11}-\sqrt{3})/2$, whereas for $\xi =(\sqrt{3}-1)/2$ it gives the SATWAP Bell inequality. In both cases the optimal CGLMP observables (expression (15) for $N=m=2$) are used.

The question we want to answer now is whether by changing ξ between the above two values we can obtain Bell inequalities maximally violated by partially entangled states (87) for various values of γ. To answer this question let us first take the observables (15) and compute the value of the Bell expression for the state (85). This gives us the following function of ξ and γ:

$$\begin{eqnarray}&&{ \mathcal J }(\xi ,\gamma )=\displaystyle \frac{4}{3}\displaystyle \frac{3+\gamma (2\sqrt{3}+\gamma -\xi \gamma )}{2+{\gamma }^{2}}.\end{eqnarray} \tag{ 88 }$$

To find its maximal value for a fixed ξ, we need to satisfy the following condition $\partial { \mathcal J }(\xi ,\gamma )/\partial \gamma =0$, which is equivalent to finding the root of a second degree polynomial in γ. The maximal value of ${ \mathcal J }(\xi ,\gamma )$ is found to be at ${\gamma }_{+}(\xi )=\left[{\left(4{\xi }^{2}+4\xi +25\right)}^{1/2}-2\xi -1\right]/2\sqrt{3}$, and it is given by

$$\begin{eqnarray}&&{{ \mathcal J }}_{\max }(\xi )=\displaystyle \frac{1}{3}\left[5-2\xi +\sqrt{25+4(\xi +1)\xi }\right].\end{eqnarray} \tag{ 89 }$$

Of course, the above derivation is not a proof that, for a given ξ, ${{ \mathcal J }}_{\max }(\xi )$ is the maximal quantum violation of the Bell inequality (85), however, based on our numerical study we conjecture this to be the case. Notice first that for $\xi =1$ and $\xi =(\sqrt{3}-1)/2$, the expression (89) reproduces the maximal quantum violations of the CGLMP and SATWAP Bell inequalities, respectively. Then, we have tested our conjecture for other values of ξ by using the NPA hierarchy, which we implemented using the Yalmip toolbox [36] and the SeDuMi solver [37] in Matlab. The NPA hierarchy provides outer approximations to the quantum set of correlations, and for a given Bell inequality, it allows one to find an upper bound on the maximal quantum violation of the Bell esxpression. We employed this technique for values of $\xi \in [-0.99,100]$ with the step 0.01, and for all these values of ξ the value obtained agrees with (89) up to solver precision $\,{10}^{-8}$, which is a strong implication that it is the maximal quantum violation of the corresponding inequality. Note that for $\xi \in [-0.99,42]$, the level $1+{AB}$ of the hierarchy was sufficient, while for $\xi \in [42,100]$ we used the level 2, except for a small amount of values in the interval $[85,100]$ for which the level $2+{AAB}$ was necessary.

In [25], we also showed how this class of inequalities can be used to self-test partially entangled states using the method of [38].

The extension of the last section to more parties turns out to be straightforward. We follow the same procedure: we start from (44) for $N=3,4$, $m=2$, $d=3$, and divide it by one the parameters so that we obtain a one-parameter class of Bell expressions

$$\begin{eqnarray}&&{J}_{N,\mathrm{2,3}}(\xi )={{\mathbb{P}}}_{0}-\xi {{\mathbb{Q}}}_{0},\end{eqnarray} \tag{ 90 }$$

with $N=3,4$ (let us notice here that both ${{\mathbb{P}}}_{0}$ and ${{\mathbb{Q}}}_{0}$ defined in equations (11) and (12) depend on N). We can compute the classical bound of these expressions, obtaining

$$\begin{eqnarray}{C}_{\mathrm{3,2,3}}(\xi )=\left\{\begin{array}{cl}-8\xi , & \mathrm{if}\,\,\xi \leqslant -1,\\ 2(3-\xi ), & \mathrm{if}\,\,-1\leqslant \xi \leqslant 1,\\ 4, & \mathrm{if}\,\,\xi \geqslant 1\end{array}\right.\end{eqnarray} \tag{ 91 }$$

and

$$\begin{eqnarray}{C}_{\mathrm{4,2,3}}(\xi )=\left\{\begin{array}{cl}-16\xi , & \mathrm{if}\,\,\xi \leqslant -10/11,\\ 10-5\xi , & \mathrm{if}\,\,-10/11\leqslant \xi \leqslant 2/5,\\ 8, & \mathrm{if}\,\,\xi \geqslant 2/5.\end{array}\right.\end{eqnarray} \tag{ 92 }$$

Let us now consider the following partially entangled GHZ states

$$\begin{eqnarray}&&| {\mathrm{GHZ}}_{\gamma }^{(N)}\rangle =\displaystyle \frac{1}{\sqrt{2+{\gamma }^{2}}}(| 0{\rangle }^{\otimes N}+\gamma | 1{\rangle }^{\otimes N}+| 2{\rangle }^{\otimes N}).\end{eqnarray} \tag{ 93 }$$

As in the previous subsection, we compute the values of ${{ \mathcal J }}^{(N)}(\xi ,\gamma )$ $(N=3,4)$ for the corresponding partially entangled GHZ states and the measurements (15) and (16), and then we solve $\partial {{ \mathcal J }}^{(N)}(\xi ,\gamma )/\partial \gamma =0$ to obtain the optimal

$$\begin{eqnarray}&&{\gamma }^{(3)}(\xi )={\gamma }^{(4)}(\xi )=\displaystyle \frac{\sqrt{4{\xi }^{2}+4\xi +25}-2\xi -1}{2\sqrt{3}},\end{eqnarray} \tag{ 94 }$$

which is the same value as for $N=2$. Substituting (94) into the values of the Bell expressions, one obtains

$$\begin{eqnarray}&&{{ \mathcal J }}_{\max }^{(3)}(\xi )=2(1+2\xi +\sqrt{25+4(\xi +1)\xi }),\end{eqnarray} \tag{ 95 }$$

and ${{ \mathcal J }}_{\max }^{(4)}(\xi )=2{{ \mathcal J }}_{\max }^{(3)}(\xi )$. We conjecture that they are the maximal quantum violations of ${J}_{\mathrm{3,3,2}}(\xi )$ and ${J}_{\mathrm{4,3,2}}(\xi )$, respectively.

To support this conjecture, we use the NPA hierarchy. With the change of scenario, it takes more time to solve each SDP, so we do not check as many values of ξ as in the section above. For $N=3$, we checked values of $\xi \in [-1,5]$ with step 0.1 and found that the gap was of order ${10}^{-7}$ or lower. For $N=4$, we checked values of $\xi \in [-1,2]$ with step 0.5 and found that the gap was of order ${10}^{-8}$ or lower.