1 Introduction↩︎

In 1964, Bell’s celebrated paper [1] transformed the long-standing debate on the completeness of quantum mechanics into a quantitative and experimentally testable problem. By deriving inequalities that must be satisfied by any theory based on local hidden variables (LHV), Bell showed that the predictions of quantum mechanics could, in principle, be distinguished from those of any locally causal model through statistical correlations observed in experiments [2].

Although originally introduced to rule out local hidden-variable descriptions of nature, Bell inequalities have gradually evolved into an essential practical resource. They now constitute the key ingredient of device-independent quantum information processing, where one aims to infer properties of quantum systems from observed measurement statistics [3]–[6]. In this paradigm, tasks such as the generation of certified randomness or the establishment of secure cryptographic keys can be guaranteed without relying on the internal workings of the quantum devices employed.

In a Bell scenario, several spatially separated parties perform measurements on subsystems of a shared physical resource and record the outcomes. Each party chooses from a set of possible measurements (inputs), and the experiment produces a joint probability distribution of outcomes conditioned on these choices. The nature of the physical theory underlying the experiment constrains the set of achievable correlations. For instance, correlations admitting a local hidden variables description must lie within a polytope defined by convex combinations of deterministic strategies [2], [7].

The hyperplanes bounding the local polytope define Bell inequalities, namely linear constraints that all LHV correlations must satisfy. While quantum mechanics allows for correlations that extend beyond this polytope and violate these inequalities, quantum correlations are themselves confined to a broader, non-polytopic convex set. The bound of a Bell inequality violation permitted by quantum theory is known as the quantum or Tsirelson bound, named after the seminal work of Boris S. Tsirelson [8]. Determining these bounds not only quantifies the gap between classical and quantum correlations but also reveals fundamental structural properties of quantum theory. A remarkable example is self-testing [4]: certain correlations residing on the border of the quantum set can essentially be performed by only one possible set of measurements and states (up to change of basis and auxiliary degrees of freedom). This is particularly relevant for certifying the correct operation of a device without the need to characterize it. Consequently, self-testing results are leveraged in many device-independent quantum security protocols across various scenarios (see e.g. [9]–[14]).

Understanding the joint geometry of local and quantum bounds is, however, non-trivial, and our current knowledge remains limited [15]. A "minimal" Bell scenario involving two users, each performing two measurements with two outputs (the \((2,2,2)\) scenario), is relatively well understood: all extremal points of the quantum correlations set have been characterized [16]–[20], and the literature offers a wide variety of Tsirelson bounds and techniques [10], [21] to obtain them using both maximally and non-maximally entangled states.

In higher-dimensional systems, as the number of measurement outcomes increases, the complexity of the quantum set scales exponentially, and analytical methods for systematically deriving sum-of-squares decompositions remain scarce. For instance, while the approach in [21] successfully derives sum-of-squares decompositions tailored to specific shared quantum states, it becomes increasingly difficult to apply as the Hilbert space dimension grows. In these complex scenarios, numerical approaches, like those based on Navascués-Pironio-Acín (NPA) hierarchy [22]–[24], are possible. However, numerical methods offer limited insight into the underlying analytical structure or the physical principles defining the quantum set. Furthermore, the NPA approach can become computationally intractable at higher levels [25] or finite levels of the hierarchy may fail to converge tightly to optimal quantum bounds in specific scenarios [26]. For these reasons, this work investigates an analytical approach to the study of the set of quantum correlations, contributing to the systematic characterization of the quantum boundary for correlations achieved using bipartite maximally entangled states in arbitrary dimensions.

From a mathematical standpoint, analytically proving Tsirelson bounds often relies on the technique of sum-of-squares (SOS) decompositions [4]. The idea is to rewrite the Bell operator \(S\) associated with a given inequality as a sum of positive semidefinite terms: \(S = \sum_i P_i^\dagger P_i\). The positivity of each term \(P_i^\dagger P_i\) immediately implies that the expectation value of \(S\) on any quantum state is lower-bounded by \(0\), thus establishing a bound. Moreover, the quantum states that achieve this bound must satisfy the constraints \(P_i \ket{\psi}=0\) for all \(i\). This can be exploited to characterize the structure of the optimal state \(\ket{\psi}\) and the corresponding measurements in self-testing [4], [19], [21], [25], [27]. In this way, the SOS method provides a powerful bridge between observed correlations and the analytical conditions required to derive robust self-testing results.

Our method builds upon these ideas and is detailed in Section 2. We begin with a "trivial" Bell operator, defined via a sum-of-squares decomposition, for which the local and quantum bounds coincide. For this baseline operator, an initial optimal strategy can be readily identified, involving a maximally entangled state and arbitrary measurements by one of the parties. The core of our approach is to then apply specific transformations to this initial operator to derive a new class of physically observable Bell operators. Crucially, we find that these transformations are a subset of unitary transformations, ensuring that the resulting Bell operator is physically observable, and constrained such that an optimal quantum strategy remains identifiable. In particular, the new optimal strategy is derived directly from the initial strategy and the applied transformation. The family of SOS Bell operators we derive have the following simple structure: \[\label{eq:SOS95intro} S = \mathbf{G}^\dagger \mathbf{G}\,,\qquad \mathbf{G}=M{\boldsymbol{\Lambda}}-U\Gamma{\boldsymbol{\Pi}} \,.\tag{1}\] This expression can be interpreted as follows. The vectors \({\boldsymbol{\Lambda}}\) and \({\boldsymbol{\Pi}}\) contain the symbolic local projectors for Alice and Bob, respectively. The matrices \(M\) and \(\Gamma\) characterize how the projectors for a given measurement enter the SOS Bell operator for Alice and Bob. Meanwhile, the unitary matrix \(U\) mixes Bob’s projectors across different measurements. Because \(S\) is a sum-of-squares decomposition, it defines a boundary \(\langle S \rangle \geq 0\). This bound is tight because it can be saturated by a maximally entangled state using a known quantum measurement strategy; thus, the condition \(\langle S \rangle \geq 0\) defines a valid Tsirelson bound. Conversely, the local bound is generally strictly greater than zero, and can be computed by standard numerical optimization or checking over the finite set of deterministic strategies, as we will discuss.

Within the compact formalism of Eq. 1 , our approach yields Bell operators that are saturated by a maximally entangled state and a generic choice of measurements (projectors) on Alice’s side, for an arbitrary number of inputs and outcomes. As special cases, we recover a large class of Bell inequalities for qubits, including those introduced in [21], [28], [29], as well as the elegant Bell inequality of [30]. For qudits, we first reproduce the Bell inequality for mutually unbiased bases (MUBs) introduced in [31] and derive a new inequality for MUBs requiring fewer inputs. We then reproduce and generalize the SATWAP inequality with two inputs per party [32].

The paper is organised as follows. In Sec. 2, we introduce the mathematical formalism and the analytical method proposed in this work. At the end of this section, we summarize the procedure in a highlighted red box. Sec. 3 is devoted to applying our method to derive quantum bounds saturated by qubit strategies. Previously known results from the literature are reproduced in highlighted blue boxes. In Sec. 4, we discuss Bell operators associated with high-dimensional, \(d >2\) strategies. We first discuss scenarios in which Alice has two inputs and Bob \(2d\) inputs, with specific examples for the MUBs, and then we consider a symmetric scenario in which both parties have two inputs. Finally, in Sec. 5, we present the outlook and discuss potential future directions of this work.

2 Description of the method↩︎

Consider a standard Bell scenario involving two spatially separated parties, Alice and Bob. Each party chooses from \(m\) possible measurement settings, denoted by \(x, y \in \{1, \dots, m\}\), which yield one of \(d\) possible outcomes, \(a, b \in \{1, \dots, d\}\). The experimental statistics are characterized by the joint probability distribution \(p(a, b | x, y)\), referred to as correlations.

In the quantum mechanical framework, we model this scenario by assigning Alice and Bob \(m\) sets of \(d\) orthogonal projectors \(\Lambda_a^{(x)}\) and \(\Pi_b^{(y)}\) acting on the Hilbert spaces \(\mathcal{H}_A\) and \(\mathcal{H}_B\), respectively. These projectors satisfy the following orthogonality and completeness relations: \[\label{eq:normalization95cons} {\Lambda}_a^{(x)} {\Lambda}_{a'}^{(x)} = \delta_{a a'} {\Lambda}_a^{(x)}, \quad \sum_{a=1}^{d} {\Lambda}_a^{(x)} = \openone, \quad \forall x\tag{2}\] and similarly for \(\Pi^{(y)}_b\). The observed correlations are obtained as the expectation value of these local operators acting on a shared quantum state \(\ket{\psi} \in \mathcal{H}_A \otimes \mathcal{H}_B\): \[\require{physics} \label{eq:correlations} p(a,b|x,y) = \expval{{\Lambda}^{(x)}_a \otimes {\Pi}^{(y)}_b}{\psi} \, .\tag{3}\]

Bell operators are observables whose expectation values are linear combinations of the correlations \(p(a, b | x, y)\). Since the range of their attainable expectation values may differ between local realistic theories and quantum theory, they can serve as witnesses that certify which physical frameworks are compatible with experimental data [2].

Identifying the theoretical boundaries of the expectation values of a Bell operator \(S_0\) is greatly simplified if \(S_0\) is written as sum-of-squares (SOS) decomposition: \[\label{eq:SOSform} S_0 = \sum_{i=1}^N P_i^\dagger P_i\,,\tag{4}\] where \(N\) is an integer and each \(P_i\) is a function of Alice’s and Bob’s local projectors, \(P_i(\Lambda_a^{(x)}, \Pi^{(y)}_b)\). The primary advantage of this decomposition is that the resulting operator is explicitly positive semi-definite. Consequently, if a specific quantum strategy yields \(\langle S_0 \rangle = 0\), then zero represents a Tsirelson bound. This quantum bound can then be compared to the classical limit, which is necessarily greater than or equal to zero.

In this work, we aim to develop a systematic method to construct Bell operators in the SOS form 4 . These operators must satisfy two requirements:

• Observability: By definition, the expectation value of a Bell operator must be expressible as a linear combination of the observed correlations 3 . This implies that \(S_0\) cannot contain products of operators from the same party with different inputs, such as \(\Lambda^{(x)}_a \Lambda^{(x')}_{a'}\) for \(x \neq x'\). Such "cross-terms" represent quantities that are not directly observable in a standard Bell scenario.

• Saturation: There must exist a quantum optimal strategy (a state \(\ket*{\widetilde{\psi}}\) and measurements \(\widetilde{{\Lambda}}^{(x)}_a\) and \(\widetilde{{\Pi}}^{(y)}_b\)) that saturates the bound \(\require{physics} \expval{S_0}=0\). This implies that \(P_i(\widetilde{\Lambda}_a^{(x)},\widetilde{\Pi}^{(y)}_b)\ket*{\widetilde{\psi}}=0\) for all \(i\).

In the remainder of this paper, we derive Bell operators whose Tsirelson bounds are saturated using maximally entangled states \(\ket*{\widetilde{\psi}}=\ket*{\Phi^+}\). Our approach involves two steps: first, we identify a simple initial Bell operator and its optimal strategy; Subsequently, we apply transformations to construct families of general operators.

A comment on notation: Throughout this paper, we denote the specific operators of a quantum strategy with a tilde (\(\tilde{\Lambda}\)), whereas the abstract symbolic operators used to define the Bell operator are written without a tilde (\(\Lambda\)).

2.1 A simple initial operator↩︎

A particularly simple way to fulfill the previous requirements is the following.

• To satisfy the observability constraint, we define each \(P_i\) as a function of only a single input \(x\) for Alice and \(y\) for Bob: \[\label{eq:no95mix95input} P_i \equiv P^{(xy)}_i = \sum_{a=1}^d p_{i,a}^{(x)}\Lambda_a^{(x)} + \sum_{b=1}^d q_{i,b}^{(y)}\Pi^{(y)}_b\tag{5}\] for some coefficients \({p_{i,a}^{(x)}, q_{i,b}^{(y)}}\). This structure ensures that the resulting SOS decomposition contains no cross-terms between different measurement settings of the same party.

• To satisfy the saturation condition, we exploit a property of maximally entangled states, \[\ket*{\Phi^+}=\frac{1}{\sqrt{d}}\sum_i \ket{i}\ket{i} \;.\] For any choice of projectors \(\widetilde{\Lambda}_a^{(x)}\) on Alice’s side, Bob’s projector \(\widetilde{\Pi}^{(x)}_a=(\widetilde{\Lambda}^{(x)}_a)^T\), where the transpose is computed in the computational basis \(\ket{i}\), satisfies (see e.g.Lemma 4.4 in [33]) \[\require{physics} \label{eq:cond95max95state} \qty(\widetilde{\Lambda}_a^{(x)} \otimes \openone-\openone \otimes \widetilde{\Pi}_a^{(x)}) \ket*{\Phi^+}=0.\tag{6}\]

  The above points tells us that if we consider the pieces of the SOS decomposition to be \[\require{physics} P_{i_x}^{(x)} = \sum_{a=1}^{d} M^{(x)}_{i_x, a} \qty(\Lambda_a^{(x)} - \Pi_a^{(x)})\,,\qquad i_x=1,\cdots, n_x\] where \(M^{(x)}\) are generic \(n_x \times d\) complex matrices, the bound \(\langle S_0 \rangle = 0\) is reached by choosing a maximally entangled state, \(\ket*{\widetilde{\psi}}=\ket*{\Phi^+}\), arbitrary sets of projectors \(\widetilde{\Lambda}_a^{(x)}\) on Alice’s side and \(\widetilde{\Pi}_a^{(x)} = (\widetilde{\Lambda}^{(x)}_a)^T\) on Bob’s. Each \(P^{(x)}_{j_x}\) represents a combination of projectors belonging to the same input \(x\), while \(n_x\) is the number of such combinations for any input \(x\).

Hence a simple Bell operator involving all the projectors and satisfying our requirements is \[\require{physics} \label{eq:basic95inequality} S_0 =\sum_{x=1}^{m}\sum_{i_x=1}^{n_x}\qty(P^{(x)}_{i_x})^\dagger P^{(x)}_{i_x}\,.\tag{7}\] To match our previous notation in Eq. 4 , we define \(N=\sum_x n_x\) and introduce a global index \(i \in \{1,\dots,N\}\). This single index replaces the nested structure, running over all different inputs and absorbing the individual sub-indices \(i_x=1,\dots,n_x\) associated with each input \(x\).

The condition \(\require{physics} \expval{S_0} = 0\) defines a Tsirelson bound, i.e., a limit on the correlations attainable within quantum mechanics. However, this bound is not particularly interesting, as it can also be saturated by local realistic strategies. In this sense, \(\require{physics} \expval{S_0} \ge 0\) constitutes both a Tsirelson inequality and a Bell inequality, making the Bell operator \(S_0\) unsuitable for probing nonlocal features of quantum mechanics or for technological applications exploiting nonlocality. Nevertheless, one may attempt to transform 7 in order to derive less trivial Bell operators. To this end, we collect the projectors into \(md\)-dimensional vectors \[\require{physics} \begin{align} {\boldsymbol{\Lambda}} &= \qty({\Lambda}^{(1)}_1 \cdots {\Lambda}^{(1)}_{d} \cdots {\Lambda}^{(m)}_1 \cdots {\Lambda}^{(m)}_{d})^T, \\ {\boldsymbol{\Pi}} &= \qty({\Pi}^{(1)}_1 \cdots {\Pi}^{(1)}_{d} \cdots {\Pi}^{(m)}_1 \cdots {\Pi}^{(m)}_{d})^T , \end{align}\] where the transpose is taken on the \(md\)-dimensional space and not on the Hilbert space. With this notation, the sum-of-squares operator 7 can be written compactly as \[\label{eq:SOS95M} S_0=({\boldsymbol{\Lambda}-{\Pi}})^\dagger M^{\dagger}M({\boldsymbol{\Lambda}-{\Pi}})\;,\tag{8}\] where we introduced a block-diagonal matrix \(M=\bigoplus_x M^{(x)}\). Note that the block diagonal structure allows to maintain the observability.

The central idea of this work is to apply suitable transformations to obtain less trivial Bell operators. Such operators will be again written in terms of an SOS decomposition, and for this reason the Tsirelson bound will be always fixed to 0, while the local bound will in general be larger than 0. Since our technique is constructed around the relation 6 , the new optimal strategies will involve the maximally entangled state. While analogous relations exist for more general entangled states, their application to constructing Bell operators is more challenging (see Sec.5).

2.1.1 SOS decompositions in terms of normal operators↩︎

Before describing our results, it is useful to discuss in more detail the expression 8 , in particular the role and meaning of the matrix \(M\). For this purpose we introduce the operators \[\begin{align} \mathbf{{A}} & \equiv M \mathbf{{\Lambda}}\,,\qquad {A}_{i_x}^{(x)} = \sum_{a}M^{(x)}_{i_x a}{\Lambda}_{a}^{(x)} \;, \\ \mathbf{{B}} & \equiv M \mathbf{{\Pi}}\,,\qquad {B}_{i_y}^{(y)} = \sum_{b}M^{(y)}_{i_y b}{\Pi}_{b}^{(y)} \;. \end{align}\] In these definitions we are combining together the projectors of a fixed input, \(x\) for Alice and \(y\) for Bob, to define sets of normal operators which satisfy the following commutation relations: \[\left[{A}_{i_x}^{(x)}, {{A}_{i_x'}^{(x)}}^\dagger\right] = \left[{B}_{i_x}^{(x)}, {{B}_{i_x'}^{(x)}}^\dagger\right] = 0 \,, \qquad \forall i_x, i'_x\] as it can be easily checked using the orthogonality relations 2 of the projectors. Being normal, for each \(x\) and \(y\) these operators can be simultaneously diagonalized for any given input and expressed as linear combinations of the original projectors. Note also that, for each \(x\) and \(i_x\), the operators \(A^{(x)}_{i_x}\) and \(B^{(x)}_{i_x}\) share the same eigenvalues. These definitions turns 8 into the diagonal form \[\begin{align} S_0&= \sum_{x=1}^m\sum_{j_x=1}^{n_x} \left( {A}_{j_x}^{(x)} - {B}_{j_x}^{(x)} \right)^\dagger \left( {A}_{j_x}^{(x)} - {B}_{j_x}^{(x)} \right)\\ &= (\mathbf{{A}}-\mathbf{{B}})^\dagger (\mathbf{{A}}-\mathbf{{B}}) =( \mathbf{{A}}^\dagger \mathbf{{A}} - \mathbf{{A}}^\dagger \mathbf{{B}} - \mathbf{{B}}^\dagger \mathbf{{A}} + \mathbf{{B}}^\dagger \mathbf{{B}} )\;. \label{eq:SOS95AB} \end{align}\tag{9}\] Using the index notation introduced above, we assign each vector \(\mathbf{A}\) (and \(\mathbf{B}\)) a global index \(j=1, \dots , N\) obtained by concatenating the \(m\) sequences of indices \(j_x\).

Parameterizing Bell operators and correlations in terms of normal operators is a well-established technique in the literature for both qubit and high-dimensional Bell inequalities [2], [4], [19], [21], [27], [29], [32], [34]. In our discussion, representing the Bell operators in terms of the set of normal operators \(\mathbf{A}\) and \(\mathbf{B}\) is a convenient way to encode the structure of the Bell operator \(S_0\) in 7 . In this sense, when we express the optimal strategy in terms of these operators \[\mathbf{\widetilde{{A}}}=M{\boldsymbol{\widetilde{{\Lambda}}}} \;, \qquad \mathbf{\widetilde{{B}}}=M{\boldsymbol{\widetilde{{\Pi}}}} \;,\] we include both the choice of initial optimal strategy and the choice of the initial operator \(S_0\).

2.2 Obtaining non-trivial Tsirelson bounds↩︎

The Bell operator \(S_0\) written in 8 defines, via its expectation value, a trivial Tsirelson bound for a scenario with \(m\) inputs and \(d\) outcomes per party, and it is obtained with the strong requirement 5 of having SOS decompositions of just one input in each term. In order to obtain more general and non-trival Bell operators from Eq.@eq:eq:SOS95M , we write the initial \(\mathbf{{A}}\) and \(\mathbf{{B}}\) sets of (symbolic) operators by means of invertible linear transformations \(V_A\) and \(V_B\) and other sets of (symbolic) operators \(\mathbf{{A}'}\) and \(\mathbf{{B}'}\): \[\label{eq:AVA} {\boldsymbol{A}} = V_A \, {\boldsymbol{A}'}, \qquad {\boldsymbol{B}} = V_B \, {\boldsymbol{B}'}\,.\tag{10}\] The role of the transformations \(V_{A/B}\), which we assume to be invertible, is to mix operators of different "initial" inputs in a controlled way. To clarify this, let’s make more explicit, in the index notations, how the primed operators 10 are defined \[\begin{align} {A}'_{k}&=\sum_{x=1}^m\sum_{i_x=1}^{n_x}(V^{-1}_{A})_{k,i_x}{{A}}_{i_x}^{(x)} \\ &= \sum_{i=1}^{N}(V^{-1}_{A})_{k,i}{\boldsymbol{A}}_i \end{align} \begin{align}\qquad \qquad {B}'_{w}&=\sum_{y=1}^m \sum_{i_y=1}^{n_y}(V^{-1}_{B})_{w,i_y}{{B}}_{i_y}^{(y)} \\ &= \sum_{i=1}^{N}(V^{-1}_{B})_{w,i}{\boldsymbol{B}}_i \;. \end{align}\] The inner sum over \(i_x\) mixes operators inside the \(x\)-th input while the external sum over \(x\) puts together contributions from different inputs. The resulting new operators \({A}'_k\) and \({B}'_w\) identify, for each \(k,w=1,\dots, N\), a new input with \(d\) outputs. So, each of them must admit a spectral decomposition in terms of its own set of projectors: \[\label{eq:diagon:Aprime} {A}'_{k}=\sum_{a=1}^d\alpha^{(k)}_{a}{{\Lambda}'}^{(k)}_{a} \;, \qquad {B}'_{w}=\sum_{b=1}^d\beta^{(w)}_{b}{{\Pi}'}_{b}^{(w)}\,.\tag{11}\] The transformed Bell operator, defined in a new scenario with \(N\) measurements per party, will be \[\label{eq:trans95sos} \begin{align} S &= {\boldsymbol{A}'}^\dagger V_A^\dagger V_A {\boldsymbol{A}'} - {\boldsymbol{A}'}^\dagger V_A^\dagger V_B {\boldsymbol{B}'} - {\boldsymbol{B}'}^\dagger V_B^\dagger V_A {\boldsymbol{A}'} + {\boldsymbol{B}'}^\dagger V_B^\dagger V_B{\boldsymbol{B}'} \, \end{align}\tag{12}\] which can be written in terms of projectors by using the relations 11 . However, we still need to satisfy our constraints of observability and saturation.

• Observability: To ensure the observability of the Bell operator \(S\), we must exclude products of operators from the same party that involve different inputs, such as \({A}'_{k}{A}'_{k'}\) for \(k \neq k'\). Looking at Eq.@eq:eq:trans95sos we see that such non-observable "cross-input" terms can only originate from the terms \(\mathbf{A'}^\dagger V_A^\dagger V_A \mathbf{A'}\) and \(\mathbf{B'}^\dagger V_B^\dagger V_B \mathbf{B'}\). When each \(k\) corresponds to a different input, in order to preserve observability, we require \(V^\dagger V\) to be diagonal for both Alice and Bob. More general expressions for \(V\) are allowed if different \(k\) can be associated with the same inputs. These aspects will be discussed in detail in the next section. Ultimately, it is sufficient to consider \(V\) to be unitary, since this choice does not limit the families of obtainable Tsirelson bounds.

• Saturation: To saturate the transformed bound \(\require{physics} \expval{S} = 0\), one might formally define an optimal strategy by considering the inverse transformations of \(V_{A/B}\) applied to the initial optimal strategies \(\mathbf{\widetilde{{A}}}\), \(\mathbf{\widetilde{{B}}}\): \[\label{eq:opti95strate95U} \mathbf{\widetilde{{A}'}} = V^{-1}_A \mathbf{\widetilde{{A}}} \;, \qquad \mathbf{\widetilde{{B}'}} = V^{-1}_B \mathbf{\widetilde{{B}}} \,.\tag{13}\] However, it is not guaranteed a priori that these transformed operators correspond to a valid quantum strategy. To be a valid quantum strategy, \(\mathbf{\widetilde{{A}'}}\) and \(\mathbf{\widetilde{{B}'}}\) must be composed of normal operators, \[\require{physics} \label{eq:commutator95condition} \qty[\widetilde{{A}'}_k, \qty(\widetilde{{A}'}_k)^\dagger] = 0, \qquad \qty[\widetilde{{B}'}_k, \qty(\widetilde{{B}'}_k)^\dagger] = 0 \quad \forall k=1, \dots,N\tag{14}\] i.e.each of their element must be diagonalizable. Therefore, to ensure the transformed Bell operator possesses a well-defined quantum bound of zero, we must identify transformations \(V_A\) and \(V_B\) that satisfy the above commutation relations. We note that, for real choices of the matrix \(M\) (i.e., when the operators \(\widetilde{{A}}_k\) and \(\widetilde{{B}}_k\) are Hermitian), a simple solution of 14 always exists, by choosing real matrices \(V_A\) and \(V_B\).

2.3 Measurement scenarios and allowed transformations↩︎

From the previous section, the Bell operator 12 is generally associated with \(N\) distinct inputs per party, one for each operator \({A}'_k\) and \({B}'_k\). However, the number of inputs can be reduced by exploiting the commutation properties of the optimal strategy that saturates the Tsirelson bound. For example, if \(\require{physics} \qty[\widetilde{{A}'}_{k_1},\widetilde{{A}'}_{k_2}]=0\), there is no need to treat \(k_1\) and \(k_2\) as separate inputs; the Tsirelson bound can be saturated by performing a single measurement defined in terms of their common projectors. Consequently, even at the level of the symbolic Bell operator, we can identify \(k_1\) and \(k_2\) as a single input. In this way, we can define a new Bell operator with fewer inputs, which effectively shares the same optimal strategy as the original operator. Both operators are equally valid, but their properties away from the saturation point (e.g., robustness to experimental noise or to finite-size effects) can generically be different, since the number of measurements is different. The optimal strategy \(\{\widetilde{{A}'}_k\}\) (\(\{\widetilde{{B}'}_k\}\)) tells us the minimal scenario to consider, namely the number \(m_A\) (\(m_B\)) of different inputs at Alice (Bob) side.

To formalize this concept, we first need to identify sets of mutually commuting operators within the final optimal strategy (the one saturating \(\ev{S}=0\)). Suppose Alice and Bob have \(m_A\) and \(m_B\) such sets, labeled by indices \(x\) and \(y\), containing \(n^A_x\) and \(n^B_y\) elements respectively. Since the operators inside each set \(x\) and \(y\) share a common set of projectors, they can be associated with the same measurement. Hence, we can adopt the notation from Sec.2.1.1, and denote an element of the sets of commuting operators as \({A}'^{(x)}_{i_{x}}\) and \({B}'^{(y)}_{i_{y}}\), where \(i_{x}=1,\dots, n^A_x\), \(i_{y}=1,\dots, n^B_y\), \(x= 1,\dots,m_A\) and \(y=1,\dots,m_B\). By applying permutations, we can always reorder \({{\boldsymbol{A}}}'\) and \({\boldsymbol{B}}'\) to group together the operators of each set \(x\) and \(y\): \[\require{physics} \begin{align} {{\boldsymbol{A}}}'&=\qty({A}'^{(1)}_{1},\cdots,{A}'^{(1)}_{n^A_1},\cdots,{A}'^{(m_A)}_{1},\cdots,{A}'^{(m_A)}_{n^A_{m_A}}) \;, \\ {{\boldsymbol{B}}}'&=\qty({B}'^{(1)}_{1},\cdots,{B}'^{(1)}_{n^B_1},\cdots,{B}'^{(m_B)}_{1},\cdots,{B}'^{(m_B)}_{n^B_{m_B}}) \;. \end{align}\] These permutations can be absorbed in \(V\) with a redefinition of the applied transformation. We can now introduce a compact expression for the spectral decomposition of commuting operators: \[\label{eq:final95spectral95decomposition} \begin{align} \widetilde{{A}'}^{(x)}_{j_{x}} & = \sum_{a=1}^d(\Gamma^{(x)}_A)_{j_{x},a}\widetilde{{\Lambda}'}^{(x)}_{a}\,, \\ \widetilde{{B}'}^{(y)}_{j_{y}} & = \sum_{b=1}^d(\Gamma^{(y)}_B)_{j_{y},b}\widetilde{{\Pi}'}^{(y)}_{b}\,, \end{align}\tag{15}\] where the rows of the matrix \(\Gamma_A^{(x)}\) represent the spectral decomposition of each operator. As we did for the matrix \(M\) before, we introduce block diagonal matrices \(\Gamma_{A/B}=\bigoplus_x \Gamma^{(x)}_{A/B}\). The spectral decomposition dictates which are the eigenvalues to be inserted in the SOS decomposition 12 to ensure a saturating quantum strategy: \[\require{physics} \label{eq:trans95sos2} \begin{align} S=\qty(V_A\Gamma_A {\boldsymbol{\Lambda}'}-V_B \Gamma_B{\boldsymbol{\Pi}'})^\dagger (V_A\Gamma_A {\boldsymbol{\Lambda}'}-V_B \Gamma_B{\boldsymbol{\Pi}'}) \,. \end{align}\tag{16}\] We emphasize that \(\Gamma\) depends strictly on the chosen measurement scenario. This \(\Gamma\) is the relevant object to consider in the presence of commutation relations in the final optimal strategy. Its dimension reflects the fact that it is not necessary to consider distinct sets of projectors if they are associated with commuting operators.

2.3.0.1 Imposing the final scenario

In principle, we could also force from the beginning the commutation relations between different primed in order to obtain a certain final Bell scenario, and this requirement restricts in general the possible \(V_{A/B}\) to be applied. For example, the requirement \(\require{physics} \qty[\widetilde{{A}'}_{k_1},\widetilde{{A}'}_{k_2}]=0\) translates into the additional constraint \[\require{physics} \label{eq:constr95commutator} \begin{align} 0 = \qty[\widetilde{{A}'}_{k_1},\widetilde{{A}'}_{k_2}] &= \sum_{s,\ell} (V_A^{-1})_{k_1s}\qty[(V_A^{-1})_{k_2\ell}]^* \,\qty[\widetilde{{\boldsymbol{A}}}_{k},\qty(\widetilde{{\boldsymbol{A}}}_{\ell})^\dagger] \\ &= \sum_{s,\ell,\mu,\nu} (V_A^{-1})_{k_1s}\qty[(V_A^{-1})_{k_2\ell}]^* \,M_{s\mu} \qty(M_{\ell,\nu})^*\,\qty[\widetilde{{\boldsymbol{\Lambda}}}_{\mu},\widetilde{{\boldsymbol{\Lambda}}}_{\nu}]\,, \end{align}\tag{17}\] for \(V_A\) to satisfy. Above condition highlights two key dependencies of the solution. First, the structure of \(V\) depends on the initial choice of Bell operator. This fact is explicit in the second line of 17 , when expressing the normal operators in terms of projectors: different choices of \(M\) may lead to distinct allowed transformations. Second, the possible solutions for \(V\) are intrinsically linked to the commutation relations of the initial optimal strategy. As we will discuss in Sec.4, certain scenarios may preclude non-trivial solutions for specific choices of \(M\) and initial strategies.

2.3.0.2 Allowed transformations \(V\)

We now consider the most general form of the transformations \(V\). For simplicity, we omit the subscripts \(A\) and \(B\), with the understanding that the following discussion applies separately to Alice and Bob. Since products of operators associated with the same inputs may appear in the Bell operator, we require that \[\label{eq:block95diag95V} V^\dagger V = \bigoplus_{x=1}^m K^{(x)} \,,\tag{18}\] where each \(K^{(x)}\) is a \(n_x \times n_x\) matrix acting between operators of the same input. By applying the polar decomposition to \(V\), we can then write: \[\require{physics} \label{eq:gen95V} V=U \Sigma \,,\qquad \Sigma=\qty(V^\dagger V)^{1/2}\,,\tag{19}\] where \(U \in \text{U}(N)\) is a unitary matrix and \(\Sigma\) Hermitian. The matrix \(\Sigma\) inherits the block-diagonal structure of \(V^\dagger V\), allowing it to be expressed in terms of a direct sum of \(m\) blocks \(\Sigma^{(x)}\), each acting exclusively on the operators of a single input: \(\Sigma = \bigoplus_x \Sigma^{(x)}\).

However, without losing generality, we could freely set \(\Sigma\) to be the identity. Indeed, since \(\Sigma\) is block-diagonal and invertible, it does not alter the commutation relations 17 solved by \(U\). Consequently, the two different transformations \(V_{1,A}=U_A\Sigma_A\) and \(V_{2,A}=U_A\) define the same final projectors \({\boldsymbol{\widetilde{{\Lambda}'}}}\). We have that \[\begin{align} {\boldsymbol{\widetilde{{A}'}}}&= V^{-1}_{1} {\boldsymbol{\widetilde{{A}}}} = \Gamma_{1,A} {\boldsymbol{\widetilde{{\Lambda}'}}}\,, \\ {\boldsymbol{\widetilde{A''}}}&= V^{-1}_{2}{\boldsymbol{\widetilde{{A}}}} = \Gamma_{2,A} {\boldsymbol{\widetilde{{\Lambda}'}}}\,, \end{align} \quad\Rightarrow\quad V_{1,A}\Gamma_{1,A}=V_{2,A}\Gamma_{2,A}\] and then \(V_{1,A}\) and \(V_{2,A}\) lead to the same final SOS 16 . Therefore, even in presence of commutation relations, the generic transformation in Eq.@eq:eq:AVA can be considered unitary without loss of generality.

2.4 Symmetries of the derivation↩︎

Before summarizing our discussion and presenting the method in a compact form, we address the symmetries of our approach. Identifying symmetries that do not alter the families of obtainable Tsirelson bounds allows us to simplify computations and restrict the set of transformations under consideration.

2.4.1 Symmetries from the choice of the initial Bell operator↩︎

The initial Bell operator \(S_0\) in 9 is defined in terms of generic normal operators, gathered in the vectors \(\mathbf{A}\) and \(\mathbf{B}\) and determined by the matrix \(M\). In general, the matrix \(M\) has no restrictions but different choices of it lead to the same family of Bell operators \(S\).

• Scale invariance: Choices of \(M\) that differ only by a global scaling factor correspond to a simple rescaling of the Bell operator \(S_0\), which is physically irrelevant. In fact, the commutation constraints are not affected, and the final optimal strategies \(\widetilde{{\boldsymbol{A}}'}\) and \(\widetilde{{\boldsymbol{B}}'}\) are invariant.

• Shift invariance: For each initial input \(x\), we can perform the shift: \[\label{eq:shift95M95id} M^{(x)} \to M^{(x)} + \mathbf{c}^{(x)} \mathbf{j}^T\tag{20}\] where \(\mathbf{c}^{(x)} = (c_1^{(x)}, \dots, c_{n_x}^{(x)})^T\) is an \(n_x\) dimensional vector and \(\mathbf{j} = (1, \dots, 1)^T\) is a \(d\)-dimensional vector. By utilizing the completeness relations \(\sum_{i_x} \Lambda_{i_x}^{(x)} = \sum_{i_x} \Pi_{i_x}^{(x)} = \openone\), we see that these transformations correspond to adding terms proportional to the identity to the corresponding normal operators: \[A_{j_x}^{(x)} \to A_{j_x}^{(x)} + c_{j_x}^{(x)}\openone, \qquad B_{j_x}^{(x)} \to B_{j_x}^{(x)} + c_{j_x}^{(x)} \openone \;.\] Since the operator \(S_0\) depends only on the differences \(A_{j_x}^{(x)} - B_{j_x}^{(x)}\), the global Bell operator remains invariant under such shifts. Moreover, since a term proportional to the identity does not affect the commutativity constraints 14 , the set of allowed unitary transformations \(U_{A/B}\) remains the same, and the transformed operators share the same projectors. Finally, for a fixed transformation \(U\), the shift 20 does not change the final Bell operators. To see this, consider the cases without and with the shift, applying the rotation only to Bob’s side for simplicity: \[\begin{align} U_B^{-1}\widetilde{{\boldsymbol{B}}}&= \Gamma_B \widetilde{{\Pi}'} \\ U_B^{-1}(\widetilde{{\boldsymbol{B}}}+\mathbf{c}\mathbf{1})&= \Gamma_B' \widetilde{{\Pi}'} \end{align} \qquad \implies \qquad \Gamma_B \widetilde{{\Pi}'}=\Gamma_B'\widetilde{{\Pi}'}-U_B^{-1}\mathbf{c}\mathbf{1}\] where we used the fact that the set of projectors of the transformed operators remains the same. The corresponding final Bell operators are: \[\begin{align} S&=(\Gamma_A \mathbf{\Lambda}'-U_B \Gamma_B \mathbf{\Pi}')^\dagger (\Gamma_A \mathbf{\Lambda}'-U_B \Gamma_B \mathbf{\Pi}') \\ S'& =(\Gamma_A \mathbf{\Lambda}'+\mathbf{c}\mathbf{1}-U_B \Gamma_B' \mathbf{\Pi}')^\dagger (\Gamma_A \mathbf{\Lambda}'+\mathbf{c}\mathbf{1}-U_B \Gamma_B' \mathbf{\Pi}')=S \;. \end{align}\] This implies that it is sufficient to consider only initial traceless normal operators.

• Unitary invariance: Since the SOS construction \(\sum_i P_i^\dagger P_i\) involves the product \((M^{(x)})^\dagger M^{(x)}\), there is the freedom to multiply \(M^{(x)}\) from the left by a unitary matrix. In other words, by applying the singular value decomposition (SVD): \[\label{eq:decomp95M} M^{(x)} = Y^{(x)} D^{(x)} W^{(x)}\,,\tag{21}\] where \(D^{(x)}\) is a positive semi-definite diagonal matrix, and \(Y^{(x)}\) and \(W^{(x)}\) are unitary matrices, we can set \(Y^{(x)} = \openone\) without loss of generality. This reduces the decomposition to: \[M^{(x)} = D^{(x)} W^{(x)} \,.\] In terms of allowed transformations, \(Y^{(x)}\) can be absorbed into the unitary matrix \(U\), which does not change the set of possible final operators.

We will further explore the utility of these symmetries in Sec.3 to simplify our computations. Their remarkable property is that they reduce the set of initial operators \(S_0\) that need to be considered, without restricting the set of final \(S\) operators that can be obtained.

2.4.2 Transforming one party is sufficient↩︎

Previously, we arrived at the expression \[S=\mathbf{G}_0^\dagger \mathbf{G}_0\,, \qquad {\boldsymbol{G}}_0=U_A \Gamma_A \boldsymbol{\Lambda}'-U_B \Gamma_B \boldsymbol{\Pi}'\,\] for the final Bell operator. Note that \(S\) remains unchanged if we apply the unitary \(U_A^\dagger\) to \({\boldsymbol{G}}_0\): \[S = {\boldsymbol{G}}_1^\dagger {\boldsymbol{G}}_1\,, \qquad {\boldsymbol{G}}_1= U_A^\dagger \mathbf{G}_0 =\Gamma_A{\boldsymbol{\Lambda}}'-U_A^\dagger U_B\Gamma_B{\boldsymbol{\Pi}}' \,.\] This expression implies that it is sufficient to consider transformations on only one of the sides. Indeed, as discussed, the matrix \(\Gamma_A\) is analogous to the matrix \(M\) described earlier. By considering the saturating strategy \(\Gamma_A \widetilde{\boldsymbol{\Lambda}'}\) for \({\boldsymbol{G}}_0\) as the initial strategy for obtaining \({\boldsymbol{G}}_1\), \(\Gamma_B\) becomes exactly the spectral decomposition that yields Bob’s operator after the transformation \(U_A^\dagger U_B\) on Bob’s side. To see this explicitly, note that the relations between the spectral decompositions of the initial operators and the transformed operators are (cf. 13 ): \[\begin{align} M \widetilde{{\boldsymbol{\Lambda}}}&=U_A\Gamma_A \widetilde{{\boldsymbol{\Lambda}}'}\,, \\ M \widetilde{{\boldsymbol{\Pi}}}&= M(\widetilde{{\boldsymbol{\Lambda}}}^T) = U_B\Gamma_B \widetilde{{\boldsymbol{\Pi}}'}\,, \end{align}\] (where the transpose acts on the Hilbert space of the operators). This also implies: \[\Gamma_B \widetilde{{\boldsymbol{\Pi}}'}=U_B^\dagger M(\widetilde{{\boldsymbol{\Lambda}}}^T)=(U_A^\dagger U_B)^{-1} \Gamma_A (\widetilde{{\boldsymbol{\Lambda}}'})^T\,.\] We can then obtain \(S\) by starting from an initial strategy with \(\Gamma_A \boldsymbol{\widetilde{\Lambda}}'\) instead of \(M \boldsymbol{\widetilde{\Lambda}}\), and applying \(U_A^\dagger U_B\) on Bob’s side. Therefore, without loss of generality, we can set \(U_A=\openone\).

2.5 Final Bell operators: summary↩︎

Summarizing the previous discussion and removing irrelevant superscripts in the final expression, the general SOS we derived is \[\label{eq:final95bell} S=\mathbf{G}^\dagger \mathbf{G}\geq0\,,\qquad \mathbf{G}=DW{\boldsymbol{\Lambda}}-U \Gamma {\boldsymbol{\Pi}}\,.\tag{22}\] In this expression we can freely choose the matrices \(D\) and \(W\) and a generic optimal strategy \(\widetilde{{\mathbf{A}}}=DW\widetilde{\boldsymbol{\Lambda}}\) on Alice’s side. The unitary matrix \(U\) must be then chosen in order to satisfy the normality constraint 14 . The optimal Bob’s strategy is obtained by applying \(U^{-1}\) to the transpose of Alice’s strategy and the matrix \(\Gamma\) is determined by the spectral decomposition of the normal operators obtained, as shown in 15 . The commutation properties of these operators dictate the number of inputs on Bob’s side.

Because Alice’s strategy can be expressed in terms of traceless operators (as shown in the previous section), Bob’s final strategy is also traceless, as it is defined as a linear combination of Alice’s operators.

2.5.1 Local bounds↩︎

While \(\require{physics} \expval{S}=0\) defines a quantum bound, the local bound, namely, the minimum expectation value achievable by local hidden variable theories, will generally be strictly greater than zero. For LHV theories, the set of correlations forms a convex polytope; thus, any local probabilistic behavior can be expressed as a convex combination of its vertices, which correspond to a finite number of deterministic strategies [2], [7]. Because the expectation value is a linear functional of the local correlations, minimizing it over the entire polytope is equivalent to performing the minimization over its vertices.

We represent the deterministic vertices using vectors of classical binary variables, \(\boldsymbol{\lambda}\) and \(\boldsymbol{\pi}\). Their elements, \(\lambda^{(x)}_a, \pi^{(y)}_b \in \{0,1\}\), define a deterministic assignment of outputs for each respective input. To ensure completeness, where exactly one outcome occurs per measurement setting, we impose the constraints: \[\sum_{a=1}^d \lambda^{(x)}_a = \sum_{b=1}^d \pi^{(y)}_b=1 \;.\] Subject to these boolean and normalization constraints, the classical variables algebraically mimic the idempotence (\(P^2 = P\)) and mutual orthogonality (\(P_i P_j = 0\) for \(i \neq j\)) of the quantum projector vectors \(\mathbf{\Lambda}\) and \(\mathbf{\Pi}\). This structural parallel allows us to compute the local bound \(S_{\rm LHV}\) directly from the sum-of-squares decomposition without expansion. The problem then reduces to minimizing the squared norm of the vector \(\mathbf{g}\) over all valid deterministic assignments: \[\label{eq:LHV95bound} \mathcal{S}_{\rm LHV} = \min_{\boldsymbol{\lambda},\boldsymbol{\pi}} \|\mathbf{g}\|^2\,,\qquad \mathbf{g} = DW\boldsymbol{\lambda}-U\Gamma \boldsymbol{\pi} \;,\tag{23}\] which gives the local bound for a given choice of initial optimal strategy and unitary transformation.

We are now ready to summarize the method discussed so far.

In the following, we apply this method to both qubit and qudit strategies. In both cases we will provide analytical solutions for the constraint in Eq.@eq:eq:recipe95constr , reproducing new results and obtaining new Bell operators with the associated Tsirelson bounds and saturating strategies.

3 Qubit inequalities↩︎

We now proceed to a detailed analysis of the Tsirelson bounds that can be saturated by generic qubit strategies acting on a maximally entangled state, in a scenario in which Alice can perform \(m\) different measurements. These strategies have two outputs for each input so that we require \(d=2\). Following the discussion of Sec.2.4, we can consider traceless normal operators for Alice’s initial optimal strategy. We note that any normal traceless \(2\times 2\) operator can be written as a phase term times a rescaled traceless unitary operator, \({A}^{(x)}_{j_x} = \exp[i \theta_{j_x}]k_{j_x}{\mathcal{A}}^{(x)}_{j_x}\). This choice aligns us with the existing literature, where qubit Bell inequalities are typically expressed in terms of unitary operators. We then write the initial Bell operator as \[S_0 = ({\boldsymbol{\Lambda}} - {\boldsymbol{\Pi}})^\dagger M^\dagger M ({\boldsymbol{\Lambda}} - {\boldsymbol{\Pi}})= ({\boldsymbol{A}} - {\boldsymbol{B}})^\dagger ({\boldsymbol{A}} - {\boldsymbol{B}}) = \frac{1}{2}({\boldsymbol{\mathcal{A}}} - {\boldsymbol{\mathcal{B}}}) K^\dagger K ({\boldsymbol{\mathcal{A}}} - {\boldsymbol{\mathcal{B}}})\,,\] where the vectors \(\boldsymbol{\mathcal{A}},\boldsymbol{\mathcal{B}}\) group the unitary operators, and the \(N\times N\) diagonal matrix \(K\) groups the coefficients \(\sqrt2k_{j_x}\). Since \(M\) is defined up to a global constant, then also \(K\) is; we therefore impose \(\require{physics} \Tr[K^2] = 1\) (\(\sum_j k_j^2 = 1/2\)).

Note that the phase terms can be removed thanks to the unitary invariance described in Sec.2.4, so Alice’s initial optimal strategy, saturating the Tsirelson bound of \(S_0\), is \[\widetilde{{\mathcal{A}}}^{(x)}_{j_x}={\vec{n} }^A_{j_x}\cdot\vec{\sigma}\,,\qquad |\vec{n}_{j_x}^A|=1\] a generic traceless unitary observable. Here we label three dimensional Euclidean vectors with an arrow, \(\vec{n}\), for differentiating from the vectors of operators defined previously. Bob’s initial optimal strategy is the transpose of Alice’s in the canonical basis \(\widetilde{{{\mathcal{B}}}}^{(x)}_{j_x}={\vec{n}}^B_{j_x}\cdot\vec{\sigma}\) with \({\vec{n}}^B_{j_x}=(n^A_{1,j_x},-n^A_{2,j_x},n^A_{3,j_x})\).

We seek the general form of an \(N \times N\) unitary matrix \(U\) satisfying the condition (cf.@eq:eq:commutator95condition ): \[\require{physics} \label{eq:qubit95m95input} \qty[\widetilde{{B}'}_j,\widetilde{{B}'}_j^\dagger]= \sum_{p, \ell = 1}^N U^{-1}_{jp} (U^{-1}_{j\ell})^* \qty[ \widetilde{{\boldsymbol{B}}}_{p}, (\widetilde{{\boldsymbol{B}}}_{\ell})^\dagger ] = \sum_{p, \ell = 1}^N U^{-1}_{jp} (U^{-1}_{j\ell})^* k_p k_\ell\qty[ \widetilde{{\boldsymbol{\mathcal{B}}}}_{p}, (\widetilde{{\boldsymbol{\mathcal{B}}}}_{\ell})^\dagger ] =0 \;.\tag{28}\] This equation requires the transformed operators, \(\widetilde{{B}'}_j = \sum_k U^{-1}_{jk} \widetilde{{\boldsymbol{B}}}_k\), to be normal. Since we started with initial traceless operators and performed linear combinations, the resulting operators will also be traceless. As already noted, any normal traceless \(2\times 2\) operator can be written as a phase term times an Hermitian traceless matrix: \(\widetilde{{B}'}_j= e^{i\theta_j}H_j = \sum_{k=1}^N U^{-1}_{jk} \widetilde{{\boldsymbol{B}}}_k\), with \(H_j = H_j^\dagger\). Multiplying both sides by \(e^{-i\theta_j}\) and requiring the right-hand side to be Hermitian, we obtain (note that \(\widetilde{{\boldsymbol{B}}}_k\) is Hermitian): \[\require{physics} \label{eq:qubit95norm95cond} \sum_{k=1}^N \qty(e^{-i \theta_j}U^{-1}_{jk} - e^{i \theta_j}(U^{-1}_{jk})^*) \widetilde{{\boldsymbol{B}}}_k = 0 \;.\tag{29}\] For every choice of initial strategy \(\widetilde{{\boldsymbol{B}}}_k\) a solution of the previous equation is found by requiring that \(e^{-i \theta_j}U^{-1}_{jk}\) are the components of a real matrix, \(R\). This is how to say that real combinations of Hermitian operators are always Hermitian. If the initial strategy consists of linearly independent operators \(\widetilde{{\boldsymbol{B}}}_k\) (e.g., the three Pauli matrices) \(e^{-i \theta_j}U^{-1}_{jk}\in \mathbb{R}\) is the only possible solution. More general cases are discussed in Appendix 6. We note that commuting normal \(2\times 2\) matrices are proportional, implying that the initial optimal strategy with more than one operator per input always has some dependencies. Such dependencies can be reduced by considering a matrix \(M\) which reduces the number of normal operators for each measurement in the initial strategy to one.

As a last step, we need to require the operator \(U^{-1}_{jk}=e^{i \theta_j}R_{jk}\) to be the components of an unitary matrix. This implies that \(U\) must be a diagonal phase matrix times an orthogonal one: \[U= \Phi O\,,\] where \(\Phi=\text{diag}(e^{-i \theta_1}, \dots , e^{-i \theta_N})\) and \(O \in {\rm SO}(N)\).

Since the matrix \(\Phi\) does not mix different inputs, as we discussed in 2.3 for the matrix \(\Sigma\), we can set \(\Phi=\openone_m\) since its contribution disappears when performing the spectral decomposition of the transformed operators. The new operators after the transformation can be written as \[\widetilde{B}'_k=\sum_{k}O^{-1}_{kj}\widetilde{\boldsymbol{B}}_j\,,\] which can be expressed in terms of unitary operators as \[\widetilde{\mathcal{B}}'_x=\frac{1}{|\vec{u}^B_x|}\vec{u}^B_x\cdot\vec{\sigma}\,,\qquad \vec{u}^B_x=\sqrt{2}\sum_{x'}O^{-1}_{xx'}\,k_{x'}\vec{n}^B_{x'}\,.\] The original Bell operator, \(S_0\), is then transformed into \[S = \frac{1}{2}\mathbf{G'}^\dagger \mathbf{G'}\,,\qquad \mathbf{G'}=K\boldsymbol{\mathcal{A}}'-O\Gamma\boldsymbol{\mathcal{B}}'\,,\] where \(\Gamma = {\rm diag}(| {\vec{u}^{B}_1}|, | {\vec{u}^{B}_2}|,\cdots, | {\vec{u}^{B}_N}|)\). By expanding the product we get \[\begin{align} S &=\frac{1}{2}(\boldsymbol{\mathcal{A}}')^\dagger K^2\boldsymbol{\mathcal{A}}'+ \frac{1}{2}\boldsymbol{\mathcal{B}}'^\dagger\Gamma^2\boldsymbol{\mathcal{B}}'- (\boldsymbol{\mathcal{A}}')^\dagger K O\Gamma\boldsymbol{\mathcal{B}}' \\ &= \openone - (\boldsymbol{\mathcal{A}}')^\dagger K O\Gamma\boldsymbol{\mathcal{B}}'\,, \end{align}\] where we computed \(\require{physics} (\boldsymbol{\mathcal{A}}')^\dagger K^2\boldsymbol{\mathcal{A}}'=\Tr[K^2] \openone = \openone\), and \(\require{physics} (\boldsymbol{\mathcal{B}}')^\dagger\Gamma^2\boldsymbol{\mathcal{B}}' = \Tr[\Gamma^2]\openone = \sum_x |\vec{u}^{B}_x|^2 \openone = \sum_x 2k^2_x|\vec{n}^{B}_x|^2 \openone = \Tr[K^2]\openone\).

3.0.0.1 Summary for qubit Tsirelson bounds

Here we summarize how our method derives Tsirelson bounds saturated by qubit strategies.

1. Choose \(N\) 3-dimensional unit vectors \(\vec{n}^A_x\), \(x=1,\cdots,N\) (\(|\vec{n}_x^A|=1\)). Define the unit vectors \({\vec{n}}^B_x=(n^A_{1,x},-n^A_{2,x},n^A_{3,x})\). These vectors correspond to \(N\) unitary operators defining Alice’s strategy. The number of sets of parallel vectors, \(m\), defines the number of inputs of Alice, each one with 2 outcomes.

2. Choose a non-negative \(N\times N\) diagonal matrix \(K\) with elements \(K_{xy}=\sqrt{2}k_x\delta_{xy}\) with the normalization condition \(\require{physics} \Tr[K^2] = \sum_x 2k_x^2=1\).

3. Choose an orthogonal matrix with unit determinant, \(O \in {\rm SO}(N)\). If the vectors \(\vec{n}^A_x\) are linearly dependent, you can also choose some non-orthogonal matrices characterized in Appendix 6.

4. Define the vectors \[\vec{u}^{B}_x=\sqrt{2}\sum_{y=1}^N (O^{-1})_{xy}\,k_y\,\vec{n}^{B}_{y} \;.\]

5. Define \({\boldsymbol{\mathcal{A}}}\) as the vector with elements \(\mathcal{A}_x\), and \(\boldsymbol{\mathcal{B}}\) the vector with elements \(\mathcal{B}_x\), where \(\mathcal{A}_x\) and \(\mathcal{B}_x\) are symbolic unitary Hermitian operators. The Bell operator is given by \[\label{eq:final95S95qubit} \begin{align} \mathbf{G'}&=K{\boldsymbol{\mathcal{A}}}-O\Gamma{\boldsymbol{\mathcal{B}}}\,, \\ S &= \frac{1}{2}\mathbf{G'}^\dagger \mathbf{G'} =\openone-{\boldsymbol{\mathcal{A}}}^\dagger KO\Gamma {\boldsymbol{\mathcal{B}}} \,,\qquad \end{align}\tag{30}\] where \[\Gamma={\rm diag}(| {\vec{u}^{B}_1}|, | {\vec{u}^{B}_2}|,\cdots, | {\vec{u}^{B}_N}|)\,.\]

6. The quantum bound \(\ev{S}\geq 0\) is saturated by the following strategy: \[\label{eq:optimal95strategy95qubit} \ket*{\Phi^+}=\frac{\ket{00}+\ket{11}}{\sqrt{2}}\,, \quad \widetilde{{\mathcal{A}}}_x=\vec{n}^A_x\cdot\vec{\sigma}\,,\qquad \widetilde{{\mathcal{B}}}_x=\frac{1}{|\vec{u}^B_x|}\vec{u}^B_x\cdot\vec{\sigma}\,.\tag{31}\]

7. The classical bound is obtained by minimizing the following function: \[\mathcal{S}_{\rm LHV}\geq \beta= \min_{\{\boldsymbol{a},\boldsymbol{b}\}} g^2({\boldsymbol{a}},{\boldsymbol{b}})\,,\qquad g^2({\boldsymbol{a}},{\boldsymbol{b}}) = 1-\,{\boldsymbol{a}}^T KO\Gamma {\boldsymbol{b}}\] where the elements \(a_x, b_x \in \{-1,+1\}\) of the vectors \(\mathbf{a}\) and \(\mathbf{b}\) define a deterministic assignment for each measurement. We note that \(g^2({\boldsymbol{a}},{\boldsymbol{b}})\) has \(2^{2N}\) discrete possible values.

3.1 Minimal scenario: two inputs per party↩︎

We now restrict our attention to the case \(N=m\), namely, the case in which there is only one normal operator per input. As noted above, this case is relevant for qubits, since commuting traceless operators (such as those associated with the same input) are essentially the same operator, up to an overall scalar factor. In the minimal scenario \((2,2,2)\), both parties have two inputs. Alice’s strategy is defined by two unit vectors \(\vec{n}^A_1\) and \(\vec{n}^A_2\) and the corresponding Bob’s initial strategy is \({\vec{n}}^B_x=(n^A_{1,x},-n^A_{2,x},n^A_{3,x})\). The two numbers \(k_{1,2}\) that satisfy \(k^2_1+k^2_2=1/2\) can be chosen as \[k_1=\frac{\cos\gamma}{\sqrt{2}}\,,\qquad k_2=\frac{\sin\gamma}{\sqrt{2}}\,,\qquad 0\leq\gamma\leq\frac{\pi}{2} \;.\] The general orthogonal matrix \(O\) can be written as \[\label{eq:rotation95O} O = \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix}\,,\tag{32}\] from which the vectors \(\vec{u}^{B}_x=\sqrt{2}\sum_{y=1}^2 (O^{-1})_{xy}\,k_y\vec{n}^{B}_{y}\) can be derived: \[\begin{align} \vec{u}^B_1&=\cos\theta \cos\gamma \,\vec{n}^B_1 - \sin\theta\sin\gamma \,\vec{n}^B_2 \;, \\ \vec{u}^B_2&=\sin\theta \cos\gamma\,\vec{n}^B_1 + \cos\theta \sin\gamma\,\vec{n}^B_2 \;. \end{align}\] The vector \({\boldsymbol{G}'}\) and the final Bell operator are then given by \[\begin{align} \mathbf{G'} &= \begin{pmatrix} \cos\gamma \mathcal{A}_1 \\ \sin\gamma \mathcal{A}_2 \end{pmatrix} - \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos \theta \end{pmatrix} \begin{pmatrix} |\vec{u}^B_1| & 0 \\ 0 & |\vec{u}^B_2| \end{pmatrix} \begin{pmatrix} \mathcal{B}_1 \\ \mathcal{B}_2 \end{pmatrix} \label{eq:G95final95trans}\,, \\ S &= \openone -\cos\theta\Bigl[\cos\gamma |\vec{u}^B_1| \mathcal{A}_1\mathcal{B}_1+\sin\gamma |\vec{u}^B_2| \mathcal{A}_2\mathcal{B}_2\Bigr]-\sin\theta \Bigl[\cos\gamma |\vec{u}^B_2| \mathcal{A}_1\mathcal{B}_2- \sin\gamma |\vec{u}^B_1| \mathcal{A}_2\mathcal{B}_1\Bigr] \;. \end{align}\tag{33}\] The Tsirelson bound \(\ev{S}\geq0\) is saturated by measuring the unitary observables \[\begin{align} \widetilde{\mathcal{A}}_1&=\vec{n}^A_1\cdot\vec{\sigma}\,,\qquad & \widetilde{\mathcal{B}}_1&=\frac{1}{|\vec{u}^B_1|} \vec{u}^B_1\cdot\vec{\sigma} \;, \\ \widetilde{\mathcal{A}}_2&=\vec{n}^A_2\cdot\vec{\sigma}\,,\qquad & \widetilde{\mathcal{B}}_2&=\frac{1}{|\vec{u}^B_2|} \vec{u}^B_2\cdot\vec{\sigma} \;. \end{align}\] The local bound can be explicitly evaluated as \[\mathcal{S}_{\rm LHV}\geq \beta =\min\{\beta_1,\beta_2\}\,,\] where \[\label{eq:local95bound95qubit} \begin{align} \beta_1&=1 -\cos\gamma \abs{\cos\theta |\vec{u}^B_1| +\sin\theta \abs{\vec{u}^B_2}} -\sin\gamma\abs{ \sin\theta |\vec{u}^B_1| - \cos\theta |\vec{u}^B_2|}\,, \\ \beta_2&=1 -\cos\gamma \abs{\cos\theta |\vec{u}^B_1| -\sin\theta \abs{\vec{u}^B_2}} -\sin\gamma\abs{ \sin\theta |\vec{u}^B_1| + \cos\theta |\vec{u}^B_2|}. \end{align}\tag{34}\]

We note that we can always choose Alice’s strategy as \(\vec{n}^A_1=(0,0,1)\) and \(\vec{n}^A_2=(\sin \alpha,0,\cos \alpha)\). This simply corresponds to fixing a global reference system for final measurement eigenvectors and the state, which remains maximally entangled.

3.2 Fixing a final optimal strategy↩︎

The previous example box also shows how to use our method to find Tsirelson bounds saturated by strategies fixed on both Alice’s and Bob’s sides. For example, suppose you want to find Bell operators whose Tsirelson bounds are saturated by the strategy \[\label{eq:max95rand95strat} \mathcal{A}_1 = \sigma_z\,,\quad \mathcal{A}_2 = \cos (a) \sigma_z + \sin (a)\sigma_x\,, \qquad \mathcal{B}_1 = \sigma_x\,,\quad \mathcal{B}_2 = \cos (b) \sigma_z - \sin (b)\sigma_x\,.\tag{37}\] This choice is relevant because the correlation between the measurements \(\mathcal{A}_1\) and \(\mathcal{B}_1\) vanishes, yielding \(\ev*{\mathcal{A}_1 \mathcal{B}_1}{\phi^+}=0\). Consequently, certifying this strategy guarantees that the measurement outcomes can be used to generate maximal global randomness in the scenario with two inputs and two outputs. A particular case is the choice \(b=\delta\) and \(a=\delta-\pi/2\), which corresponds to the measurements saturating 36 and are self-tested in case where \(\ev{I_\delta}=\mathcal{I}_{\mathcal{Q}}'\).

To check the possibility of deriving a boundary saturated by the previous measurements, we solve the system of equations \[(K{\boldsymbol{\mathcal{A}}})_j = (O \Gamma {\boldsymbol{\mathcal{B}}})_j^T \quad \forall j\,,\] in terms of \(K,O\) and \(\Gamma\). A solution is found by setting \[\require{physics} \begin{align} \gamma &= \operatorname{arccot}\qty(-\frac{\cos(a)\sin(a+b)}{\sin(b)})\,, \\ \theta &= \arctan(-1-\cot(b)\tan(a))\,, \end{align} \qquad \begin{align} |\vec{u}_1^B| &= \tan(b)\sqrt{-\cot(b)\tan(a+b)}\,,\\ |\vec{u}_2^B| &= \sqrt{\frac{\cos(a)\cos(a+b)}{\cos(b)}}\,, \end{align}\] provided that \(-\pi/2 < a+b<0\). The parameters defined above characterize a family of Tsirelson bounds saturated by the strategy in 37 , and may therefore be utilized in device-independent protocols for global randomness generation; a formal study of their certification properties remains an objective for future works.

3.3 Asymmetric scenarios and generalization of the elegant Bell inequality↩︎

Here we leverage the discussion of Sec.2.3 to retrieve the Bell operator of the Elegant Bell Inequality (EBI) [30], defined in an asymmetric scenario, as a particular case of a more general and much broader family of operators with a symmetric number of inputs for Alice and Bob.

We consider the case in which Alice’s measurements consist of four qubit observables forming the vertices of a regular tetrahedron on the Bloch sphere \[\label{eq:tetrahedron} \begin{align} \widetilde{{\mathcal{A}}}_1 &= \frac{\sigma_x - \sigma_y + \sigma_z}{\sqrt{3}}\equiv \vec{n}_1^A \cdot \vec{\sigma}\\ \widetilde{{\mathcal{A}}}_2 &= \frac{\sigma_x + \sigma_y - \sigma_z}{\sqrt{3}}\equiv \vec{n}_2^A \cdot \vec{\sigma} \end{align}\qquad \begin{align} \widetilde{{\mathcal{A}}}_3 &= \frac{-\sigma_x - \sigma_y - \sigma_z}{\sqrt{3}} \equiv \vec{n}_3^A \cdot \vec{\sigma}\\ \widetilde{{\mathcal{A}}}_4 &= \frac{-\sigma_x + \sigma_y + \sigma_z}{\sqrt{3}} \equiv \vec{n}_4^A \cdot \vec{\sigma} \end{align}\tag{38}\] satisfying the linear dependence: \(\sum_{k=1}^4 \widetilde{{\mathcal{A}}}_k = 0\). These measurements are known to be included in the optimal strategy of the EBI [35].

We take the above measurements as initial ones for our procedure, we define \(\widetilde{{\mathcal{B}}}_k \equiv (\widetilde{{\mathcal{A}}}_k)^T\) and we consider the family of Bell operators defined by four dimensional orthogonal transformations, \(O\). As discussed in Sec.3, orthogonal transformations are always valid, also in the case in which the initial measurements are linearly dependent. When we consider the \(K\) contribution and we apply \(O^{-1}\), we get the new operators \[\require{physics} \widetilde{{B}'}_j = \sum_{\ell=1}^4 (O^{-1})_{j\ell}\, \sqrt{2} k_\ell\widetilde{{\mathcal{B}}}_\ell = \qty[\sum_{\ell=1}^4 (O^{-1})_{j\ell}\, \sqrt{2} k_\ell\,\vec{n}^B_\ell] \cdot \vec{\sigma} \equiv \vec{u}_j \cdot \vec{\sigma}\,\] with \(\vec{u}_j\) a real not normalized vector. The Bell operator we find is then written as \[S = \openone - \sqrt{2} \sum_{\mu,\ell=1}^{4} k_\mu O_{\mu,\ell}|\vec{u}_\ell|\,\mathcal{A}_\mu \mathcal{B}_\ell\,.\] This Bell operator is defined in a symmetric scenario where Alice and Bob both have four dichotomic measurements with two possible outcomes \(\pm1\). The Bell operator is saturated \(\langle S\rangle=0\) when Alice chooses the measurements given in 38 , while Bob chooses the measurement \(\widetilde{{\mathcal{B'}}}_\ell=\frac{1}{|\vec{u}_\ell|}\vec{u}_\ell\cdot\vec{\sigma}\).

To modify the above operator for an asymmetric scenario in the number of inputs between Alice and Bob, we need to impose the commutator \[\require{physics} \qty[\widetilde{{B}'}_j, (\widetilde{{B}'}_k)^\dagger]=\qty[\vec{u}_j \cdot \vec{\sigma}, \vec{u}_k \cdot \vec{\sigma}] = 2i(\vec{u}_j \times \vec{u}_k)\cdot \vec{\sigma} =0\] for some choice of \(j,k\). As discussed in Sec.2.3, the number of measurements is determined by the number of sets of commuting operators in the final optimal strategy. From this we conclude that \(O\) is characterized by the equation \(\vec{u}_j \times \vec{u}_k=0\), which implies that the two vectors are collinear or that at least one of them is the zero vector. A simple choice reducing Bob’s number of measurements is to satisfy directly the linear dependence described below 38 and choose \(O^{-1}\) with one row proportional to the vector \((1, 1, 1, 1)\). This would produce a null operator in the final optimal strategy. More general solutions also exist, and one such example gives precisely the EBI, presented in the box at the end of this section.

We finally note that we can at most find a pair of operators to commute. Indeed, the transformation to the new vectors \(\vec{u}_{j}\) is governed by a \(4 \times 4\) orthogonal matrix. Because \(O\) is full-rank (it is invertible), it preserves the dimension of the space. Therefore, the four new vectors \(\{\vec{u}_{1}, \vec{u}_{2}, \vec{u}_{3}, \vec{u}_{4}\}\) must also span a 3-dimensional space, as the original \(\vec{a}_k\). If we had two separate pairs of commuting operators (e.g., \(\vec{u}_1 \propto \vec{u}_2\) and \(\vec{u}_3 \propto \vec{u}_4\)), all four vectors would be constrained to just two lines, which contradicts the requirement that the vectors must span a 3-dimensional space.

4 Higher-dimensional inequalities↩︎

As the dimension \(d\) of the Hilbert space for the optimal strategy increases, the landscape of available quantum correlations becomes significantly more complex, as do the potential relations between operators of the same or different inputs. Hence, complete and fully general arguments become difficult to carry out. Nevertheless, in this section, we discuss solutions that are valid for any number of Alice’s inputs and any choice of her strategy. We then restrict our focus to the case of two inputs (\(m=2\)), while also introducing constraints on Bob’s side.

4.1 General solutions↩︎

Let us first recap what we found in the qubit case. In \(d=2\), any traceless normal operator is proportional to an Hermitian traceless operator. By using the unitary invariance, we can eliminate irrelevant phases and choose Alice’s strategy as made exclusively with Hermitian operators. When we want to construct nontrivial Bell operators, we combine the transpose of these initial Hermitian operators (which are still Hermitian). Now, a generic linear combination of independent Hermitian operators remains Hermitian if and only if the coefficients are real. In the bidimensional case, this requirement restricted the transformations to real coefficients (up to potential linear dependencies) when combining operators from different inputs in the initial optimal strategy to identify a new input in the final optimal strategy. For this reason, to obtain relevant Bell operators, we restricted the transformations from the unitary group to the real orthogonal group.

In \(d\) dimensions, the structure is more complex. Each traceless normal operator can be expressed as a complex combination of \((d-1)\) commuting, traceless Hermitian operators. This can be seen as a direct consequence of the spectral decomposition and the completeness relation of the set of projectors, and allows in principle to a much richer structure than the \(d=2\) case. Simple arguments like the one around equation 29 do not apply, and in principle also complex transformations play a role.

However, while the \(d\)-dimensional structure allows for more intricate algebraic relations, it remains true that a real linear combination of Hermitian operators is itself Hermitian (and therefore diagonalizable), representing a well-defined measurement input. Consequently, we can generalize the qubit construction and find a solution in \(d\) dimensions considering real matrices \[(DW)_{ij} \in \mathbb{R} \quad \text{and} \quad U\equiv O \in \text{O}(md)\,,\] to obtain the valid sum of squares decomposition: \[\label{eq:qudit95gen} S = \mathbf{G}^\dagger \mathbf{G}\,, \qquad \mathbf{G}=DW{\boldsymbol{\Lambda}}-O\Gamma{\boldsymbol{\Pi}}\,,\tag{39}\] whose Tsirelson bound is saturated by a generic choice of Alice’s strategy \(\widetilde{{\mathbf{\Lambda}}}\) and a corresponding strategy for Bob \(\widetilde{{\mathbf{\Pi}}}\) defined by the projectors of the operators \(\widetilde{{\boldsymbol{B}}}=O^{-1} (\mathbf{\widetilde{{\Lambda}}}^T)\) (where the transpose is on the operators acting on the Hilbert space of the state). This formula is the high-dimensional analogue of the recipe we gave for the qubit case in Sec. 3. In that case, we also used the Bloch vector representation to simplify our results: the generalization to \(d\) dimensions requires substituting the Pauli matrices with the \(d^2-1\) generalized Gell-Mann matrices. It is important to note, however, that unlike the qubit case, not every vector in \(\mathbb{R}^{d^2-1}\) corresponds to an operator with a valid physical spectrum, as the geometry of the state space becomes more complex in higher dimensions [36].

4.1.0.1 More scenarios

The number of inputs in Bob’s scenario depends on the choice of matrix \(DW\) and on the commutation relations of the final optimal strategy, as we already discussed in Sec. 2.3. This number can be reduced if the transformed operators commute: \[\require{physics} \qty[\widetilde{{B}'}_k, \widetilde{{B}'}_{k'}] = 0.\] In \(d=2\) this can only happen if the two traceless operators are proportional, while in \(d>2\) two traceless normal operators can commute even without being proportional. However, there are some constraints we need to consider.

First, the set of normal operators that define the initial strategy spans a space in the set of matrices acting on a \(d\)-dimensional space. When we apply an invertible transformation to that set, the dimension of the spanned space cannot change.² This means, for example, that we cannot hope to have proportional operators in the transformed strategy unless there are some linear dependencies in the initial one. Second, for some choices of initial strategies, it is not even possible to construct two or more non-proportional commuting normal operators. We show this in the particular case of some mutually unbiased bases in Appendix 7. The combination of these aspects restrict the possibility of obtaining particular scenarios on Bob’s side, when Alice’s strategy is fixed.

In the following we discuss some particular cases of the general Bell operator expressed in 39 . We start by showing, in the following box, that the inequality introduced in [31] is included in this general formalism. To the best of our knowledge, the Bell inequalities of [31] are the only ones in the literature which are maximally violated by two MUBs in arbitrary dimension \(d>2\) [37]. In the next section, we will propose alternatives with fewer inputs.

4.2 Two inputs on Alice’s side↩︎

This section provides an extended treatment of the scenario with two inputs on Alice’s side. We first consider a simple application of formula 39 , leading to a scenario with \(2d\) inputs on Bob’s side. We discuss in particular detail the case in which Alice’s strategy is constructed from MUBs providing a Bell operator with fewer inputs than the one of [31]. Then, we discuss more general solutions than 39 that reduce the number of Bob’s inputs, allowing us to reproduce and generalize the SATWAP operator for two inputs.

4.2.1 \(m_A=2\) and \(m_B=2d\)↩︎

We consider formula 39 choosing \(DW=\openone\), so that Alice’s optimal strategy simply consists of two sets (\(m_A=2\)) of \(d\) orthogonal rank-1 projectors, with \(d >2\). We apply a simple \(2d\)-dimensional rotation:\[O = \begin{pmatrix}\cos\theta \openone_d & -\sin\theta \openone_d \\ \sin\theta \openone_d & \cos\theta \openone_d\end{pmatrix}\] and obtain the SOS decomposition \[\label{eq:S95rot95d} S = \sum_{k=1}^{d} \left[ \left( \Lambda_k^{(1)} - \cos\theta B_k + \sin\theta B_{k+d} \right)^2 + \left( \Lambda_k^{(2)} - \sin\theta B_k - \cos\theta B_{k+d} \right)^2 \right]\tag{43}\] where each \(B_k\) or \(B_{k+d}\), with \(k=1, \dots, d\), is an operator associated with an input. Expanding this SOS we find: \[2\openone-\sum_{k=1}^{d} \Bigg[-B_k^2 - B_{k+d}^2 + 2\left(\cos\theta \Lambda_k^{(1)} + \sin\theta \Lambda_k^{(2)}\right)B_k + 2\left(-\sin\theta \Lambda_k^{(1)} + \cos\theta \Lambda_k^{(2)}\right)B_{k+d} \Bigg] \succeq 0\;.\] The Bell inequality is completely defined when the eigenvalues of the operators \(B_k\) are specified. As discussed in the general method, every choice of Alice’s optimal projectors \(\boldsymbol{\widetilde{{\Lambda}}}\) corresponds to a determined choice of the \(B_k\)’s eigenvalues, as detailed below. Denoting Bob’s initial strategy by \(\widetilde{{\boldsymbol{\Pi}}}\), which is the transpose of Alice’s, his final optimal operators are given by: \[\label{eq:transf95op95qudit} \begin{align} \widetilde{{B}'}_k &= \cos\theta \, \widetilde{{\Pi}}^{(1)}_k + \sin\theta \, \widetilde{{\Pi}}^{(2)}_k \\ \widetilde{{B}'}_{k+d} &=-\sin\theta \, \widetilde{{\Pi}}^{(1)}_{k} + \cos\theta \, \widetilde{{\Pi}}^{(2)}_{k} \end{align}\,, \quad\qquad k=1,\cdots, d \;.\tag{44}\] Because these operators are constructed from two rank-1 projectors, their total rank can never exceed 2, regardless of the overall dimension \(d\). This means that each normal operator has only \(2\) non-vanishing eigenvalues. By defining \[\require{physics} \Tr(\widetilde{{\Pi}}^{(1)}_j \widetilde{{\Pi}}^{(2)}_j) = \Tr(\widetilde{{\Lambda}}^{(1)}_j \widetilde{{\Lambda}}^{(2)}_j) = \abs{c_{j}}^2\] the eigenvalues of the operators \(\widetilde{{B}'}_k\) and \(\widetilde{{B}'}_{k+d}\) are found to be \[\require{physics} \label{eq:eigen95qudit} \begin{align} \lambda^{(k)}_{\pm}&=\frac{1}{2} \left( \cos\theta + \sin\theta \pm \sqrt{1 - \sin(2\theta)\qty(1 - 2\abs{c_k}^2)} \right) \\ \lambda^{(k+d)}_{\pm}&=\frac{1}{2} \left( \cos\theta - \sin\theta \pm \sqrt{1 + \sin(2\theta)\qty(1 - 2\abs{c_{k}}^2)} \right) \end{align}\tag{45}\] where \(\lambda^{(k)}_{\pm}\) (\(\lambda^{(k+d)}_{\pm}\)) are the two non-vanishing eigenvalues of \(\widetilde{{B}'}_k\) (\(\widetilde{{B}'}_{k+d}\)) and \(k=1,\cdots,d\). If we now express every operator \(B_k\) in its spectral decomposition \(B_k=\sum_b\lambda^{(k)}_b\Pi^{(k)}_b\), the above Bell inequality can be written in terms of the symbolic projectors \(\{\Lambda^{(x)}_a,\Pi^{(k)}_b\}\): \[\label{eq:qudit95expanded} \begin{align} S=2\openone-2&\sum_{k=1}^{d}\sum_{b=\pm} \Big[\lambda^{(k)}_b \left(\cos\theta \Lambda_k^{(1)} + \sin\theta \Lambda_k^{(2)}-\frac{1}{2}\lambda^{(k)}_b \right)\Pi^{(k)}_b + \\ & +\lambda^{(k+d)}_b \left(-\sin\theta \Lambda_k^{(1)} + \cos\theta \Lambda_k^{(2)}-\frac{1}{2}\lambda^{(k+d)}_b \right)\Pi^{(k+d)}_b \Big] \succeq 0\;. \end{align}\tag{46}\]

4.2.1.1 The MUB case

A simple application of 46 is obtained by choosing two mutually unbiased bases on Alice’s side. The defining geometric property of MUBs in \(d\) dimensions is that the overlap probability scales inversely with the dimension, so \(\abs{c_k}^2= 1/d\). With this substitution, the eigenvalues of the two blocks in 45 become independent of the input \(k\): \[\require{physics} \begin{align} \lambda_{\pm}^{(k)}(\theta)&= \frac{1}{2} \left( \cos\theta + \sin\theta \pm \sqrt{1 - \sin(2\theta)\qty(1 - 2/d)} \right) \;, \qquad \\ \lambda_{\pm}^{(k+d)}(\theta)&= \frac{1}{2} \left( \cos\theta - \sin\theta \pm \sqrt{1 + \sin(2\theta)\qty(1 - 2/d)} \right)\,. \end{align}\]

Note that our construction yields a Bell operator with two inputs and \(d\) outcomes on Alice’s side, and \(2d\) inputs and \(d\) outcomes on Bob’s side (of which \(d-2\) are null outcomes). This result can be compared with the one in [31] (reproduced in equation 40 ), which instead requires \(d^2\) inputs on one of the two sides. The (im)possibility of constructing Bell operators with even fewer inputs on Bob’s side for generic MUBs is discussed in Appendix 7.

For a better comparison with the results in the literature, we also need to understand the relation between quantum and classical bounds of 46 . For simplicity we only discuss the case \(\theta=\pi/4\). Define \[\label{eq:T} I=2\openone-S \;.\tag{47}\] From the SOS decomposition, we know that the quantum bound of \(I\) is \(\mathcal{I}_\mathcal{Q} = 2\) while the local bound is given by (see Appendix 8) \[\mathcal{I}_{\rm LHV} = \max \left[ \frac{3}{2} + \frac{1}{\sqrt{d}} - \frac{1}{2d}, 2\sqrt{1 - \frac{1}{d}} \right]\,,\] where the first term is larger when \(d\le 5\) and the second is larger when \(d \ge 6\). Since both terms are smaller than \(\mathcal{I}_\mathcal{Q}=2\), this result shows the possibility, in principle, to use the Bell operator 43 for proper quantum protocols.

Figure 1 displays the quantum-to-local bound ratio, \(\mathcal{I}_{\mathcal{Q}}/\mathcal{I}_{\rm LHV}\), for the Bell operators defined in Eqs. 41 (from Ref. [31]) and 47 . We find that, for all \(d \ge 3\), the operator in Eq. 47 achieves a strictly larger ratio.

4.2.2 The \(m_A=m_B=2\) symmetric scenario↩︎

We now try to find Bell operators in which Bob’s strategy also involves two inputs. This investigation is driven by two primary motivations. First, this specific configuration has been previously explored in the literature (e.g., [27], [32]). We aim to verify if these results can be reproduced and extended with our techniques (a question we answer in the affirmative). Second, this scenario yields transformation examples that go beyond the orthogonal ones previously introduced and it allows us to examine the requirement of commutativity between sets of transformed variables. Beyond the specific findings, this serves to show some techniques that can be used to derive explicit transformation matrices when specific constraints are placed on Bob’s scenario.

Since this section is technically dense, its structure is summarized here to facilitate navigation through the results. The analysis begins by identifying the specific measurements on Alice’s side required to obtain the symmetric scenario, leading to Equation 48 and to the choice of the CGLMP measurements in Eq.@eq:eq:Alice95CGLMP . After reformulating the problem into a more convenient framework in Eqs.@eq:eq:cons951 and 59 , the generalized version of the SATWAP operator for two inputs is derived in Eq.@eq:eq:gen95SATWAP95op .

Comment on the notation: In this last section we will not need to use different subscripts for normal operators and for projectors. We will instead use letters \(i,j,k, \dots\) in the range from \(0\) to \(d-1\) (instead of \(1\) to \(d\)) to match the existing literature and simplifying the notation.

4.2.2.1 Choosing Alice’s measurements

As we learn from the previous discussion and Appendix 7, a generic strategy on Alice’s side will in general lead to a Bell operator with \(2d\) inputs on Bob’s side without the possibility of reducing this number to only two inputs. So, we need to look for some guidance to understand which specific strategy to choose for Alice to achieve this scenario.

Remember that the operators in Bob’s optimal strategy are combinations of the transpose of the \(2d\) projectors, \(\widetilde{{\Lambda}}_j^{(x)}=\ketbra*{\alpha^{(x)}_j}\), which define Alice’s strategy. To have only two inputs these operators must form two sets of \(d\) commuting operators and a natural idea is to require, for each set, these operators to be independent, which is typical for self-testing strategies. This means that they can be combined to produce a rank-one projector \(\widetilde{{\Pi}'}=\ketbra*{v}\) that can be written as \[\widetilde{{\Pi}'}=\ketbra*{v}=\sum_j b_j^{(1)} \ketbra*{\alpha_j^{(1)}}+\sum_j b_j^{(2)} \ketbra*{\alpha_j^{(2)}}\] for some coefficients \(b_j^{(x)}\). Now, we project this operator equation into the mixed basis by sandwiching it between \(\bra*{\alpha_k^{(1)}}\) and \(\ket*{\alpha_{\ell}^{(2)}}\): \[\braket*{\alpha_k^{(1)}}{v}\braket*{v }{\alpha_\ell^{(2)}} = b_k^{(1)} \braket*{\alpha_k^{(1)}}{\alpha_\ell^{(2)}} + b_\ell^{(2)} \braket*{\alpha_k^{(1)}}{\alpha_\ell^{(2)}} \;.\] Factoring out the overlap matrix element \(c_{k\ell}^{(12)} \equiv \braket*{\alpha_k^{(1)}}{\alpha_\ell^{(2)}}\), we find: \[\label{eq:cauchy95form} c_{k\ell}^{(12)} = \frac{\braket*{\alpha_k^{(1)}}{v}\braket*{v}{\alpha_\ell^{(2)}}}{b_k^{(1)} + b_\ell^{(2)}} \;.\tag{48}\] The expression we just found has the structure of a Cauchy-like matrix: the dependence on the rows and columns factorizes in the numerator, while in the denominator it appears through a sum. So, to obtain the scenario we have in mind, with two inputs on both parties, and \(d\) independent operators for each input of Bob’s optimal strategy, it is necessary that the two bases associated with Alice’s optimal strategy are related by a Cauchy-like transformation. Moreover, \(c_{ij}\) must also be unitary, being a change of basis between two orthonormal bases. A remarkable example satisfying these requirements is given by the measurements used in the CGLMP inequalities [38]: \[\label{eq:Alice95CGLMP} {\widetilde{{A}}_{j}^{(x)}}=\sum_{k=0}^{d-1}\omega^{kj} \ketbra*{\alpha^{(x)}_k} \;, \qquad \ket*{\alpha^{(x)}_k}=\frac{1}{\sqrt{d}}\sum_{q=0}^{d-1}\omega^{(k-\frac{x}{2}+\frac{1}{4})q}\ket*{q}\tag{49}\] where \(\omega=e^{2\pi i/d}\) and \(x=1,2\), which is characterized by the change of basis matrix \[c_{jk}^{(12)} =\frac{1}{d} \sum_{q=0}^{d-1} \omega^{-(j-k + 1/2)q}=\frac{2}{d}\left(\frac{1}{1-\omega^{k-j-1/2}}\right)\] with the manifest Cauchy-like structure. Considering this as Alice’s strategy, we now discuss which class of Bell operators can be found.

4.2.2.2 Rephrasing the commutation problem

As said, we want two sets of \(d\) commuting operators and, up to an irrelevant permutation, we can always assume that the commutation happens between the first \(d\) and the second \(d\) operators. This structure mirrors the commutation property of Bob’s initial strategy, \(\widetilde{{B}}_j^{(x)}=(\widetilde{{A}}_j^{(x)})^T\), and motivates a small change in our notation. We will divide the unitary matrix \(U\) into four blocks and equip it with two kinds of indices so that we will write it \((U)_{ij}^{(xy)}\). Indices on the top \(x,y=1, \dots m\) are related to the inputs and indices on the bottom \(i,j=0,\dots,d-1\) are associated with different outputs of the same input. With this notation the transformed operators in 10 are \[\widetilde{{B}'}^{(y)}_k=\sum_{y'}\sum_{k'}(U)_{kk'}^{(yy')}{\widetilde{{B}}}_{k'}^{(y')}\] and the commutation constraints are \[\require{physics} \label{eq:comm95eq95for95SATWVAP} \qty[\widetilde{{B}'}_k^{(y)},\qty(\widetilde{{B}'}_{k'}^{(y)})^\dagger]= \sum_{i,j=0}^{d-1} \sum_{x, z=1}^2 U_{ki}^{(yx)}\qty(U_{k'j}^{(y z)})^* \qty[\widetilde{{B}}_i^{(x)},\qty(\widetilde{{B}}_j^{(z)})^\dagger] = 0 \;,\tag{50}\] where the \(\widetilde{{B}}^{(x)}_j\) are the transposes of the matrices \(\widetilde{{A}}_j^{(x)}\) in 49 . Note that the matrix \(M\), defining the operators \(\widetilde{{B}}_j^{(x)}\) in terms of projectors, is (up to an irrelevant rescaling) unitary: \[M^{(x)}_{jk}=\omega^{jk} \;.\] Hence it can be just interpreted as a change of basis from the projectors to the normal operators. Since we can formulate the problem in the more convenient basis, we choose the projector basis and perform the change of basis at the end. So, in equation 50 , we will replace the normal operators \(\widetilde{{B}}^{(x)}_i\) and \(\widetilde{{B}}^{(z)}_j\) with the corresponding projectors.

The central objects to analyze are the commutators between the projectors of the two bases of the initial strategy. Setting \(\widetilde{{\Pi}}_i^{(x)}=\ketbra*{\beta_i^{(x)}}\): \[\require{physics} \qty[\widetilde{{\Pi}}_i^{(x)},\widetilde{{\Pi}}_j^{(z)}]=c_{ji}^{(zx)}\ketbra*{\beta_i^{(x)}}{\beta_j^{(z)}}-c_{ij}^{(xz)}\ketbra*{\beta_j^{(z)}}{\beta_i^{(x)}} \;,\] where we note that, since \(\widetilde{{\Pi}}_i^{(x)}=(\widetilde{{\Lambda}}_i^{(x)})^T\), then \(\braket*{\beta_i^{(x)}}{\beta_j^{(z)}}=\braket*{\alpha_j^{(z)}}{\alpha_i^{(x)}}={c_{ji}^{(zx)}}\). Considering the bracket with \(\bra*{\beta_{\alpha}^{(x)}}\) and \(\ket*{\beta_{\beta}^{(x)}}\) (\(\alpha,\beta=0,\dots,d-1\)), gives the components \[\require{physics} \label{eq:comm95comp} \bra*{\beta_{\alpha}^{(x)}}\qty[\widetilde{{\Pi}}_i^{(x)},\widetilde{{\Pi}}_j^{(z)}]\ket*{\beta_{\beta}^{(x)}} = c_{ji}^{(zx)}c_{\beta j}^{(xz)} \delta_{\alpha i} - c_{ij}^{(xz)} c_{j \alpha}^{(zx)} \delta_{i \beta}\equiv (C_{\alpha \beta})_{ij}^{(xz)} \;.\tag{51}\] So we can reformulate the problem as \[\require{physics} \label{eq:constr95commutator95m2} \sum_{i,j=1}^d \sum_{x, z=1}^2 U_{ki}^{(yx)}\qty(U_{k'j}^{(y z)})^* \qty(C_{\alpha \beta})_{ij}^{(xz)} = 0 \;.\tag{52}\] In matrix form, this choice allows us to write the constraint 52 as the requirement that the transformed matrices have vanishing block-diagonal components: \[\require{physics} \label{eq:comm95const95gen} \qty( U \, (C_{\alpha \beta}) \, U^\dagger )^{(yy)}_{kk'} = 0, \qquad \forall\, y,k,k',\alpha,\beta.\tag{53}\] Using the algebra of \(U(dm)\), the transformation can also be expressed in terms of a Hermitian matrix \(X\) such that \(U=e^{i X}\): \[\require{physics} \label{eq:const95exp} \qty(e^{iX} \, (C_{\alpha\beta}) \, e^{-iX})^{(yy)}_{kk'} =0 \qquad \forall\, y,k,k',\alpha,\beta.\tag{54}\] We decompose the matrix \(X\) into a block-diagonal and a block anti-diagonal component, \[X=\begin{pmatrix} X^{(11)} & X^{(12)} \\ X^{(21)} & X^{(22)} \end{pmatrix} = \begin{pmatrix} X^{(11)} & 0 \\ 0 & X^{(22)} \end{pmatrix} + \begin{pmatrix} 0 & X^{(12)} \\ X^{(21)} & 0 \end{pmatrix} \equiv X_1+X_2 .\] Here \(X^{(11)}\) and \(X^{(22)}\) are Hermitian matrices, while \(X^{(12)}=(X^{(21)})^\dagger\) is arbitrary. Each block \(X^{(xy)}\) carries indices \(i,j\). To simplify the problem, we make the ansatz that the block-diagonal and anti-block-diagonal parts commute, \[\label{eq:ansatz95comm} [X_1,X_2]=0 \;.\tag{55}\] Under this assumption, \[e^{iX} \, (C_{\alpha \beta}) \, e^{-iX} = e^{iX_1}e^{iX_2} \,(C_{\alpha \beta}) \, e^{-iX_2}e^{-iX_1}.\] Since \(X_1\) does not affect the block structure, its action can be ignored for our purposes and can in fact be reabsorbed as we discussed for the matrix \(\Sigma\) in 19 . We therefore focus on the transformation generated by \(X_2\), which we expand as \[\require{physics} \label{eq:series} e^{iX_2} \,(C_{\alpha \beta}) \, e^{-iX_2} = \sum_{n=0}^{\infty} \frac{\qty[iX_2,(C_{\alpha \beta})]_n}{n!},\tag{56}\] where \([X,Y]_n = \overbrace{[X,\ldots,[X,[X,Y]]]}^{n\;\text{times}}\) denotes the \(n\)-fold nested commutator and \([X,Y]_0 = Y\). Writing explicitly the term \(n=1\) of the series (\(n=0\) is automatically solved), we obtain \[\label{eq:first95cond} [X_2,(C_{\alpha \beta})]= \begin{pmatrix} X^{(12)}(C_{\alpha \beta})^{(21)}-(C_{\alpha \beta})^{(12)}X^{(21)} & 0 \\ 0 & X^{(21)}(C_{\alpha \beta})^{(12)}-(C_{\alpha \beta})^{(21)}X^{(12)} \end{pmatrix}.\tag{57}\] Imposing that the diagonal blocks vanish forces the entire commutator to vanish. In this case, all higher-order nested commutators in 56 also vanish, ensuring that the full series has vanishing block-diagonal components and thus satisfies the constraint. Therefore, it suffices to solve the following set of linear equations: \[\begin{align} X^{(12)}(C_{\alpha \beta})^{(21)}-(C_{\alpha \beta})^{(12)}X^{(21)}&=0, \tag{58}\\ X^{(21)}(C_{\alpha \beta})^{(12)}-(C_{\alpha \beta})^{(21)}X^{(12)}&=0,\tag{59} \end{align}\] for every \(\alpha\) and \(\beta\), which run over the elements of the operator basis used to decompose the commutators of the initial strategy (and in principle the choice of basis can be different for the two sets of equations). These equations impose constraints on the coefficients of the matrix \(X_2\). Once a solution for \(X_2\) is obtained, we can obtain \(U\) through exponentiation. Since the matrix \(X_2\) is block anti-diagonal, we have a closed expression for its exponential, in terms of cosines and sines.

4.2.2.3 Allowed transformations

By calling \(x_{ij}\) the variables of \(X^{(12)}\) and using 51 , the constraints 58 and 59 become \[\begin{align} &(x_{i\alpha}^*-x_{i \beta}^*)c_{j \alpha}^*c_{j \beta}+(x_{j \alpha}-x_{j \beta})c_{i \alpha}^*c_{i \beta}=0 \;, \\ &(x_{\alpha j}-x_{ \beta j})c_{\alpha i}^*c_{ \beta j}+(x_{\alpha i}^*-x_{ \beta i}^*)c_{\alpha j}^*c_{ \beta j}=0 \;, \end{align}\] where \(c^{(12)}_{jk}\equiv c_{jk}\). The first immediate results are obtained by choosing \(i=j\), which gives \[\label{eq:constraint95real} \Re(x_{i \alpha})=\Re(x_{i \beta})\,, \;\qquad \Re(x_{\alpha i})=\Re(x_{\beta i})\,,\tag{60}\] namely the real part of \(x_{ij}\) is constant in the entire matrix, and it is a free parameter. Hence, the previous equations can be written in terms of the imaginary parts, let us use the letter \(y\) to denote them: \[\begin{align} &(y_{i\alpha}-y_{i \beta})c_{j \alpha}^*c_{j \beta}-(y_{j \alpha}-y_{j \beta})c_{i \alpha}^*c_{i \beta}=0 \;, \\ &(y_{\alpha j}-y_{ \beta j})c_{\alpha i}^*c_{ \beta j}-(y_{\alpha i}-y_{ \beta i})c_{\alpha j}^*c_{ \beta j}=0 \;. \end{align}\] Focus on the first equation, and consider a \(j\) for which the coefficients \(c_{j\alpha}^*c_{j\beta}\) are non zero. We find \[y_{i\alpha}-y_{i \beta}=\frac{y_{j \alpha}-y_{j \beta}}{c_{j\alpha}^*c_{j\beta}}c_{i \alpha}^*c_{i \beta} \;.\] In order to satisfy this equation, the fraction must be independent of \(j\), and equal to a parameter function only of \(\alpha\) and \(\beta\), call it \(\lambda_{\alpha \beta}\). Repeating the same procedure also for the second equation, we finally get: \[\label{eq:constraint95imaginary95text} y_{i \alpha}-y_{i \beta}=\lambda_{\alpha\beta}c_{i \alpha}^*{c_{i \beta}} \;, \qquad y_{\alpha i}-y_{\beta i}=\lambda'_{\alpha\beta}c_{\alpha i}^*c_{\beta i}\tag{61}\] for some parameters \(\lambda_{\alpha \beta}\) and \(\lambda'_{\alpha \beta}\), which is analogous to the requirement we found for the real part.

4.2.2.4 Bell operators

Conditions 61 can be explicitly solved (see Appendix 9). It is simpler to represent the result by performing a permutation of the basis, and moving from the one in which the initial vector is made of \((B_0^{(1)}, \dots, B_{d-1}^{(1)},B_0^{(2)}, \dots, B_{d-1}^{(2)})\) to \((B_0^{(1)}, B_{0}^{(2)},\dots, B_{d-1}^{(1)}, B_{d-1}^{(2)})\). In this basis, the unitary matrix is a direct sum: \[U=\bigoplus_{k=1}^{d-1} \begin{pmatrix} \cos \beta & i e^{-i \pi k / d} \sin\beta \\ ie^{i \pi k / d} \sin \beta & \cos\beta \end{pmatrix}\] where we assume the interval \(0\le \beta \le \pi/2\). When \(\beta=0\) or \(\beta=\pi/2\) the operators collapse to a combination of operators of the first or second input, respectively. The transformed Bell operator is \[\require{physics} \label{eq:gen95SATWAP95op} \begin{align} S=\sum_{k=1}^{d-1}&\qty[{A}_k^{(1)}-e^{i\beta \left( 1 - \frac{2k}{d} \right)} \qty(\cos \beta {B}_k^{(1)}-i e^{-i\pi k/d}\sin \beta {B}_k^{(2)})]^2 +\\ &+\qty[{A}_k^{(2)}-e^{i\beta \left( 1 - \frac{2k}{d} \right)} \qty(-i e^{i\pi k/d}\sin \beta {B}_k^{(1)}+\cos \beta {B}_k^{(2)})]^2 \;, \end{align}\tag{62}\] where \({A}_k^{{(x)}}\) and \({B}_k^{{(y)}}\) are unitary operators with eigenvalues \(\omega^{jk}\), with \(j=0,\dots, d-1\). The bound \(\require{physics} \expval{S}=0\) is saturated by choosing \(\widetilde{{A}}_j^{(x)}\) as in 49 and \[\begin{align} {\widetilde{{B}}_k}^{(1)} &= \sum_{j=0}^{d-1} \omega^{jk} \ketbra{u_j(\beta)}{u_j(\beta)} \\ \ket{u_j(\beta)} &= \frac{1}{\sqrt{d}} \sum_{q=0}^{d-1} \omega^{-\left( j - \frac{1}{4} - \frac{\beta}{\pi} \right) q} \ket{q} \end{align} \;, \qquad \begin{align} {\widetilde{{B}}_k}^{(2)} & = \sum_{j=0}^{d-1} \omega^{jk} \ketbra{v_j(\beta)}{v_j(\beta)} \\ \ket{v_j(\beta)} &= \frac{1}{\sqrt{d}} \sum_{q=0}^{d-1} \omega^{-\left( j - \frac{3}{4} - \frac{\beta}{\pi} \right) q } \ket{q} \end{align} \;.\] By expanding the SOS and multiplying by \(1/2\) we obtain \[\label{eq:gen95satwap95exp} \begin{align} I\equiv\sum_{k=1}^{d-1} e^{-i\beta\left(1 - \frac{2k}{d}\right)} \Big[& \cos\beta {A}_k^{(1)} {B}_{d-k}^{(1)} + i e^{i\pi k/d} \sin\beta {A}_k^{(1)} {B}_{d-k}^{(2)} \\ &+ i e^{-i\pi k/d} \sin\beta {A}_k^{(2)} {B}_{d-k}^{(1)} + \cos\beta {A}_k^{(2)} {B}_{d-k}^{(2)} \Big]\,, \end{align}\tag{63}\] where we used \((B^{(y)}_k)^\dagger=B^{(y)}_{d-k}\). The quantum bound of \(I\) is \(\mathcal{I}_Q=2(d-1)\openone\), while the local bound is \[\require{physics} \mathcal{I}_{\rm LHV} = \begin{cases} \frac{\sin(2\beta)}{2} \qty[ 2\cot\left(\frac{\beta}{d}\right) - \cot\left(\frac{\beta + \pi/2}{d}\right) + \cot\left(\frac{\pi/2 - \beta}{d}\right) ] - 2 & \text{for } 0 \le \beta \le \pi/4 \\ \\ \frac{\sin(2\beta)}{2} \qty[ \cot\left(\frac{\beta}{d}\right) + 2\cot\left(\frac{\pi/2 - \beta}{d}\right) - \cot\left(\frac{\pi - \beta}{d}\right) ] - 2 & \text{for } \pi/4 < \beta \le \pi/2 \end{cases} \;.\] The family of inequalities 63 generalizes the SATWAP inequality, recovered in the special case \(\beta=\pi/4\). In Fig. 2 we show the ratio between the quantum and local bounds of the Bell operator 63 as a function of \(\beta\) and for several values of \(d\).

5 Conclusions and outlook↩︎

In this work, we addressed the problem of determining Tsirelson and local bounds for arbitrarily high-dimensional systems, in the case of two separated parties. We presented a systematic derivation of Bell operators whose Tsirelson bounds are saturated by bipartite symmetric strategies and maximally entangled states in arbitrary dimensions. Specifically, we represent Bell parameters \(S\) in terms of quadratic forms and, by applying suitable transformations, derive Tsirelson and local bounds for which an optimal strategy is guaranteed to exist. This existence is ensured by the fact that the transformations map sets of normal commuting operators to other sets of normal commuting operators. Furthermore, these transformations ensure the SOS remains observable within a Bell scenario, implying that the resulting quantum bounds are suitable for device-independent protocols, such as DI quantum key distribution or DI quantum random number generation. As examples, in our discussions we derived several known results, both for qubit and qudit strategies, and extended them with new families of inequalities.

This work provides a foundation upon which several promising research directions and open questions can be explored. First, an important open question concerns the completeness of this method in characterizing the quantum boundary. Although the recovery of several known and novel results suggests considerable generality, establishing the formal limits of our transformation-based approach remains an important challenge. Moreover, we have not yet investigated sufficient or necessary conditions on Alice’s measurements and unitary transformations \(U\) under which the Tsirelson bounds differ from the local bounds.

Beyond these considerations, future work may proceed along the following directions:

• Deriving SOS saturated by optimal strategies using non maximally entangled states.

• Extending the framework to systems involving more than two parties (like, for example, those considered in [39]–[41]).

• Systematically investigating which bounds possess robust self-testing properties.

• Identifying specific device-independent quantum protocols using our results.

We are planning to discuss some aspects of multipartite scenarios and self-testing properties in following works. Here we briefly discuss some ideas and problems that arise for non-maximally entangled states.

5.0.0.1 Non-maximally entangled states

Our current construction begins with a trivial Bell operator, constructed such that a maximally entangled state \(\ket*{\phi^+}\) naturally resides in its kernel. This construction relied on the fundamental relation for maximally entangled states: \[\label{eq:trick95max95ent} (A \otimes \openone)\ket{\phi^+} = (\openone \otimes A^T)\ket{\phi^+}\tag{64}\] which allows us to saturate a bound of \(0\) using normal operators on both Alice’s and Bob’s sides. For a full-rank non-maximally entangled state \(\ket{\psi} = \sum_i \lambda_i \ket{ii}\) an analogous relation exists. Define the diagonal matrix of Schmidt coefficients \(\Lambda = \sum_i \lambda_i \ketbra{i}{i}\). Then for every Alice operator \(A\), there exists a Bob operator \(B\) (again the transpose of Alice’s) such that \[(A \otimes \openone - \openone \otimes \Lambda A^T \Lambda^{-1})\ket{\psi} = 0 \;.\] While this provides a mathematical identity to "move" an operator from Alice to Bob, it introduces some problems. The effective operator \(B_{ef}=\Lambda A^T \Lambda^{-1}\) appearing on Bob’s side is generally not normal for non-maximally entangled states. Consequently, we cannot simply construct a physical Bell operator in the form \(S = (A - B_{ef})^\dagger (A - B_{ef})\) where both \(A\) and \(B_{ef}\) represent standard local measurements. The non-normality of \(B_{ef}\) implies that it does not correspond to a single physical observable in the Bell scenario. For this reason, all the discussion presented in this paper cannot be directly applied. Moreover, it could be the case that the sum-of-squares method is not the optimal choice for deriving Tsirelson bounds for non-maximally entangled states [25]. Still, the central question and idea remain: can we construct general, physically valid Bell operators and Tsirelson bounds for non-maximally entangled states by applying unitary or more general transformations to an initial trivial Bell operator?

Acknowledgments↩︎

LC thanks Flavio Baccari and Tommaso Grigoletto for useful discussions. This work was supported by European Union’s Horizon Europe research and innovation program under the project Quantum Secure Networks Partnership (QSNP), grant agreement No 101114043. Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or European Commission-EU. Neither the European Union nor the granting authority can be held responsible for them.

6 Linearly dependent qubit strategies↩︎

In this appendix, we investigate the effect of linear dependencies among Alice’s operators on the set of transformations that yield non-trivial Bell operators. We concentrate on the qubit case, for which the constraints on the unitary transformations \(U\) were derived in the main text: \[\require{physics} \sum_{k=1}^m \qty(e^{-i \theta_j}U_{jk} - e^{i \theta_j}U_{jk}^*) \widetilde{{B}}_k = 0.\] If Alice’s operators are linearly independent (e.g., the three Pauli matrices) the same holds for \(\widetilde{{B}}_k\), and the only way the sum can equal zero is if every individual coefficient is exactly zero. This means that \(e^{-i \theta_j}U_{jk}\) must be the components of a real matrix. In general, however, the initial set of \(m\) matrices may have \(n\) independent linear dependencies, so that there are \(n\) different vectors \(v^{(1)}, v^{(2)}, \dots, v^{(n)}\) such that, for any of them, \[\label{eq:lin95dep} \sum_{k=1}^m v_k^{(\ell)} \widetilde{{B}}_k = 0 \;, \quad \text{for } \ell \in \{1, \dots, n\} \;.\tag{65}\] This allows for a non-trivial imaginary part: \[\require{physics} U_{jk} = e^{i\theta_j}\left(R_{jk} + i \,\text{Im}\qty[\sum_{\ell=1}^n \eta_{j\ell} v_k^{(\ell)}]\right) \;,\] where \(R_{jk}\) is the purely real part, and \(\eta_{j \ell}\) are completely arbitrary parameters. This last degree of freedom acts like a sort of "gauge" term, in the sense that it does not affect the final expression of Bob’s optimal strategy, thanks to 65 . However, they may enter in the actual expression of the Bell operator .

When requiring \(U\) to be unitary, if \(\eta_{j\ell}=0\) (which is always a valid solution) then (ignoring the phases) \(U\) must be orthogonal, as discussed in the main text. If, instead, there are linear dependencies and it’s possible to choose \(\eta_{j\ell} \neq 0\), the requirement of unitarity is more involved. In this case, the matrix \(U\) can be written compactly as: \[U = R + iN\] up to a irrelevant phase-diagonal matrix and where \(\require{physics} N_{jk} = \text{Im}\qty[\sum_{\ell=1}^n \eta_{j\ell} v_k^{(\ell)}]\). For \(U\) to be a unitary matrix, it must satisfy \(U^\dagger U = \openone\): \[U^\dagger U = (R^T - iN^T)(R + iN)= \openone \;.\] Since both \(R\) and \(N\) are strictly real matrices, we find the conditions \[\label{eq:qubit95lin95indip} R^T R + N^T N = \openone \;, \qquad R^T N - N^T R = 0 \;.\tag{66}\] So, when there is linear dependence between the initial operators, the general solution also admits a non-orthogonal \(R\), with a non-trivial relationship between the real part \(R\) and the arbitrary parameters in \(\eta\) and \(v\).

7 Restrictions for MUBs↩︎

The number of inputs on Bob’s side in 43 could, in principle, be reduced using the commutation properties of the operators in 44 , but this cannot be freely done for a generic pair of MUBs, as we now show. Suppose Alice’s initial strategy consists of the measurement projectors associated with two mutually unbiased bases. We will restrict our analysis to a specific standard choice. To do so, we introduce the generalized Pauli operators \(X\) and \(Z\) acting on a local \(d\)-dimensional Hilbert space, defined by their action on the computational basis \(\{|0\rangle, |1\rangle, \dots, |d-1\rangle\}\): \[X\ket*{j} = \ket*{(j+1) \pmod d}, \qquad Z\ket*{j} = \omega^j\ket*{j}\] where \(\omega = e^{2\pi i / d}\). These operators satisfy the Weyl commutation relation: \[Z^k X^m = \omega^{km} X^m Z^k \;.\] For the MUBs in Alice’s strategy, we choose the eigenbases of \(Z\) and \(X\), and restrict our focus to the case where the dimension \(d\) is a prime number.

Suppose we want to construct two operators, in our notation they would be \(\widetilde{{\Pi}'}_1\) and \(\widetilde{{\Pi}'}_2\), which are linear combinations of the projectors of these two MUBs. We can take them to be traceless. We are saying that they can be written as \[\widetilde{{\Pi}'}_1=Z_1+X_1 \qquad \widetilde{{\Pi}'}_2=Z_2+X_2\] where \(Z_j\) are some combinations of the projectors of the first basis, and \(X_j\) of the second basis. Because the \(d^2\) matrices \(X^m Z^k\) form a complete orthogonal basis for the space of operators, any traceless operator diagonal in the first basis can be written purely as powers of \(Z\), and any traceless operator diagonal in the second basis can be written purely as powers of \(X\): \[Z_1 = \sum_{k=1}^{d-1} a_k Z^k \;, \qquad X_1 = \sum_{m=1}^{d-1} b_m X^m \;, \qquad Z_2 = \sum_{k=1}^{d-1} c_k Z^k \;, \qquad X_2 = \sum_{m=1}^{d-1} d_m X^m \;.\]

We require the mixed operators to commute: \([Z_1 + X_1, Z_2 + X_2] = 0\), which simplifies to \([Z_1, X_2] = [Z_2, X_1]\). Let’s compute the commutator \([Z_1, X_2]\): \[[Z_1, X_2] = \sum_{k,m} a_k d_m [Z^k, X^m] = \sum_{k,m} a_k d_m (\omega^{km} - 1) X^m Z^k \;.\] Evaluating \([Z_2, X_1]\) in the same way and equating the two expressions, we obtain: \[\sum_{k,m} a_k d_m (\omega^{km} - 1) X^m Z^k = \sum_{k,m} c_k b_m (\omega^{km} - 1) X^m Z^k \;.\] Because the operators \(X^m Z^k\) are linearly independent basis vectors, their coefficients must match exactly term-by-term. Therefore, for every \(k\) and \(m\): \[a_k d_m (\omega^{km} - 1) = c_k b_m (\omega^{km} - 1) \;.\] Since \(d\) is prime, and \(k, m\) are strictly between \(1\) and \(d-1\), the product \(km\) is never a multiple of \(d\). Therefore, \(\omega^{km}\) is never equal to \(1\). Because \((\omega^{km} - 1) \neq 0\), we can safely divide both sides by it, and we are left with: \[a_k d_m = c_k b_m \;.\] If one of the two final operators is a combination of only one set of projectors, say \(d_m=0\), then this implies that the other operator must also be a combination of the same set (\(b_m=0\)). Suppose instead \(d_m \neq 0\). Then we divide by it: \[a_k = \left(\frac{b_m}{d_m}\right) c_k \;.\] Let \(\lambda = b_m / d_m\). We now have \(a_k = \lambda c_k\) for all \(k\), which dictates \(Z_1 = \lambda Z_2\). Similarly, \(X_1 = \lambda X_2\). Therefore, in \(d\) dimensions with \(d\) prime, linear combinations of the projectors of the MUBs associated with \(Z\) and \(X\) commute if and only if they are proportional.

These findings restrict the scenarios that can be considered when looking for sum-of-squares decompositions where one party’s strategy consists of a pair of generic mutually unbiased bases in arbitrary dimensions [31]. Globally, the two sets of projectors of two generic MUBs span a \(2d-1\) dimensional space. Applying the unitary transformation yields \(2d\) normal operators, which define \(2d\) distinct inputs on Bob’s side. To reduce the number of these inputs and simplify the scenario, some of these operators must commute. For this to hold for a generic pair of MUBs in an arbitrary dimension, the commuting operators must be proportional. However, no more than two operators can be proportional, as exceeding this limit would reduce the dimension of the generated space too much.

8 A \(d\)-dimensional local bound↩︎

In this appendix we compute the local bound of the operator \[2\openone - S=\sum_{k=1}^{d} \Bigg[-B_k^2 - B_{k+d}^2 + \sqrt{2}\left( \Lambda_k^{(1)} + \Lambda_k^{(2)}\right)B_k +\sqrt{2}\left(\Lambda_k^{(1)} - \Lambda_k^{(2)}\right)B_{k+d} \Bigg]\] with \(d >2\), introduced in Sec.4.2. Here, \(\Lambda_k^{(x)}\) are two sets of projectors and Bob’s operators have \(d-2\) zero eigenvalues while the others are: \[\label{eq:lambda95bob} \lambda_{\pm}= \begin{cases} \lambda_{\pm}^{(1)}\equiv\frac{1}{\sqrt{2}} \left(1\pm \frac{1}{\sqrt{d}} \right) \qquad &1 \le k \le d \\ \lambda_{\pm}^{(2)}\equiv \pm \frac{1}{\sqrt{2}} \sqrt{1- \frac{1}{d}} \qquad &d+1 \le k \le 2d \end{cases}\tag{67}\]

In the local scenario, we treat all operators as commuting scalar variables. For deterministic classical strategies, Alice’s measurement outputs exactly one deterministic value for each input. This means replacing the operator \(\Lambda_k^{(x)}\) with a boolean variable \(\lambda_k^{(x)}\in \{0, 1\}\), and for a fixed \(x\), exactly one \(k\) evaluates to \(1\) while the rest are \(0\). Bob’s outcomes instead are encoded in the variable \(b_k\) and \(b_{k+d}\), restricted to the eigenvalues 67 and the \(d-2\) null eigenvalues. We want to maximize the scalar sum: \[\require{physics} I_{\rm LHV} = \sum_{k=1}^d \left[ -b_k^2 + \sqrt{2}\qty(\lambda_k^{(1)}+\lambda_k^{(2)}) b_k - b_{k+d}^2 + \sqrt{2}\qty(\lambda_k^{(1)}-\lambda_k^{(2)}) b_{k+d} \right] \;.\] Because \(\sum_k \lambda_k^{(x)} = 1\), Alice has two distinct deterministic behaviors to consider over the inputs:

• Case A: Alice outputs the same \(k\) for both input \(1\) and input \(2\). For that specific \(k\), \(\lambda_k^{(1)} = \lambda_k^{(2)} = 1\). To maximize Bob’s side, we assign \(b_{k+d} = 0\) and pick the largest available eigenvalue for \(b_k\), which is \(\lambda_+ = \frac{1}{\sqrt{2}}(1 + \frac{1}{\sqrt{d}})\). Substituting these in gives the maximum for this scenario: \[\mathcal{I}_{\rm LHV}= -\lambda_+^2 + 2\sqrt{2}\lambda_+ = \frac{3}{2} + \frac{1}{\sqrt{d}} - \frac{1}{2d} \;.\]

• Case B: Alice outputs \(k_1\) for input \(1\) and \(k_2\) for input \(2\) (\(k_1 \neq k_2\)). For the term \(k_1\) in the sum, we pick \(b_{k_1} = \lambda_+^{(1)}\) and \(b_{k_1+d} = +\frac{1}{\sqrt{2}}\sqrt{1 - 1/d}\). The maximum for this term is \(\sqrt{1 - 1/d}\). For \(k_2\), the maximum similarly evaluates to \(+\sqrt{1 - 1/d}\). Summing both contributions yields the maximum for this case: \[\mathcal{I}_{\rm LHV} = 2\sqrt{1 - \frac{1}{d}} \;.\]

The local bound is whichever strategy yields the largest number: \[\mathcal{I}_{\rm LHV} = \max \left[ \frac{3}{2} + \frac{1}{\sqrt{d}} - \frac{1}{2d}, 2\sqrt{1 - \frac{1}{d}} \right] \;.\] Which one is larger depends on the dimension \(d\): For lower dimensions (\(d \le 5\)), the first strategy is the maximum. For higher dimensions (\(d \ge 6\)), the second one.

9 Generalization SATWAP↩︎

In this appendix, as a case study, we want to show the details of the derivation and generalization of SATWAP results [32]. Here, we present all the details not explicitly written in Sec.4.2.2.

From the discussion in the main text, we need to solve, in \(y\), the equations: \[\label{eq:constraint95imaginary} y_{i \alpha}-y_{i \beta}=\lambda_{\alpha\beta}c_{i \alpha}^*{c_{i \beta}} \;, \qquad y_{\alpha i}-y_{\beta i}=\lambda'_{\alpha\beta}c_{\alpha i}^*c_{\beta i}\tag{68}\] for some parameters \(\lambda_{\alpha \beta}\) and \(\lambda'_{\alpha \beta}\) and where \[c_{jk}=\frac{2}{d}\left(\frac{1}{1-\omega^{k-j-1/2}}\right) \;.\] The basis transition matrix satisfies some properties which significantly simplify the problem, namely \[\begin{align} \label{eq:overlap95cond} c_{i \alpha}^*{c_{i \beta}}=\xi_{\alpha \beta}\left(h_{\alpha i}-h_{\beta i} \right) \\ c_{\alpha i}^*{c_{\beta i}}=\xi'_{\alpha \beta}\left(h'_{\alpha i}-h'_{\beta i} \right) \end{align}\tag{69}\] with \[\begin{align} \xi_{\alpha \beta} &= \frac{4}{d^2} \frac{1}{1 - \omega^{\beta-\alpha}} \\ h_{\alpha j} &= \frac{1}{1 - \omega^{j-\alpha+1/2}} \end{align} \;, \qquad \qquad \begin{align} \xi'_{\alpha \beta} &= \frac{4}{d^2} \frac{1}{1 - \omega^{\alpha-\beta}} \\ h'_{\alpha j}& = \frac{1}{1 - \omega^{\alpha-j+1/2}}=h_{j\alpha} \end{align} \;.\] Let us now focus on the left equation in 68 . Using the factorization property the (trivial) identity \[(y_{i \alpha}-y_{i \beta })+(y_{i \beta}-y_{i \gamma })+(y_{i \gamma}-y_{i \alpha})=0\] can be written as \[\label{eq:lin95ind} h_{\alpha i}(\lambda_{\alpha \beta} \xi_{\alpha \beta} - \lambda_{\gamma \alpha} \xi_{\gamma \alpha})+h_{\beta i}(\lambda_{\beta \gamma} \xi_{\beta \gamma}-\lambda_{\alpha \beta} \xi_{\alpha \beta}) + h_{\gamma i}(\lambda_{\gamma \alpha} \xi_{\gamma \alpha}-\lambda_{\beta \gamma} \xi_{\beta \gamma})= 0 \;.\tag{70}\] Since the \(h_{\alpha i}\) form a linearly independent basis, we obtain that the coefficients in equation 70 must vanish, namely that \(\lambda_{\alpha \beta}\xi_{\alpha \beta}=C\), for a constant \(C\). Substituting this back, \[y_{i j}=y_{i 0}+\frac{C}{\xi_{j0}}c_{i j}^*{c_{i 0}}=y_{i0}+C(h_{ji}-h_{0i}) \;.\] The first column, \(y_{i0}\), can be arbitrarily chosen (\(d\) variables), together with the (real) scaling constant \(C\). The second constraint in 68 can be tackled in exactly the same way, leading to \[y_{j i}=y_{0 i}+\frac{C'}{\xi'_{j0}}c_{ j i}^*{c_{0 i}}=y_{0i}+C'(h'_{ji}-h'_{0i})=y_{0i}+C'(h_{ij}-h_{i0})\] with a constant \(C'\), and this time we can arbitrarily choose the first row, \(y_{0i}\). However, we now need to simultaneously solve these two sets, so that \(C\) and \(C'\) are not independently free parameters. To find a closed form for any \(y_{ij}\), we first express the boundary elements \(y_{i0}\) and \(y_{0j}\) in terms of the single scalar \(y_{00}\). By choosing \(i=0\) in both equations, we have: \[\begin{align} y_{0j} &= y_{00} + C(h_{j0}-h_{00}) \;, \\ y_{j0} &= y_{00} + C'(h_{0j}-h_{00}) \;. \end{align}\] We can now plug back, obtaining two closed expressions \[\begin{align} y_{ij} &= y_{00}+C'(h_{0i}-h_{00})+C(h_{ji}-h_{0i}) \;, \\ y_{ij} &= y_{00}+C(h_{j0}-h_{00})+C'(h_{ji}-h_{j0}) \;, \end{align}\] which imply the consistency condition: \[C' \left( h_{0i}-h_{00}-h_{ji}+h_{j0}\right) = C \left( h_{j0}-h_{00}-h_{ji}+h_{0i} \right) \qquad \Rightarrow \qquad C=C' \;.\] Therefore, the general elements of the matrix are \[\label{eq:imag95const} y_{ij} = y_{00}+C(h_{ji}-h_{00}) \;.\tag{71}\] Note that the quantity \((h_{ji}-h_{00})\) is pure imaginary, so also \(C\) must be imaginary to satisfy 71 .

The \(y\)’s we just derived are the imaginary parts of the elements of the matrix \(X^{(12)}\), introduced in the main text. We already argued, around Eq.@eq:eq:constraint95real , that the real parts of these elements are constants. Combining them to form the full complex matrix elements, and rescaling \(i C \to b\), with \(b \in \mathbb{R}\), we find \[x_{jk} = a+b h_{kj}\] where \(a=x_{00}+iy_{00}-iCh_{00}\) is a generic complex constant and \(h_{kj}=1/(1-\omega^{j-k+1/2})\). Note that the constant term corresponds to a shift in the direction of the identity.

Change of basis↩︎

As we said, the variables \(x_{jk}\) are elements of a generic matrix \(X^{(12)}\), and we need to consider the exponential of \[\label{eq:X295app} X_2= \begin{pmatrix} 0 & X^{(12)} \\ {X^{(12)}}^\dagger & 0 \end{pmatrix} \;.\tag{72}\] Being block off-diagonal, the exponential can be computed. However we can make the problem even easier by changing basis and moving to the basis of the normal operators used to saturate the SATWAP, by computing \(M X_2 M^{-1}\) with \(M=\bigoplus_xM^{(x)}\) and \[M^{(x)}_{jk} = \omega^{jk} \;.\] The transformation of each block is given by \(P = M^{(x)} X^{(12)} (M^{(x)})^{-1}\). Using that \((M^{(x)})^{-1}_{kn} = \frac{1}{d} \omega^{-kn}\), the elements of the matrix \(P\) are: \[p_{mn} = \sum_{j=0}^{d-1} \sum_{k=0}^{d-1} M^{(x)}_{mj} x_{jk} (M^{(x)})^{-1}_{kn} = \frac{1}{d} \sum_{j=0}^{d-1} \sum_{k=0}^{d-1} \omega^{mj} \left( a + \frac{b}{1 - \omega^{j-k+1/2}} \right) \omega^{-kn} \;.\] We can separate this into a constant term \(p^{(a)}_{mn}\) and a shifted term \(p^{(b)}_{mn}\), corresponding to the two elements in the brackets. For the constant term \(a\), we use the orthogonality of the roots of unity, \(\sum_{k=0}^{d-1} \omega^{ck} = d \delta_{c,0}\): \[p^{(a)}_{mn}= \frac{a}{d} \left( \sum_{j=0}^{d-1} \omega^{mj} \right) \left( \sum_{k=0}^{d-1} \omega^{-kn} \right)= a d \delta_{m,0} \delta_{n,0} \;.\] For the term involving \(b\), we substitute \(s = j-k\) (so \(j = k+s\)). Taking the indices modulo \(d\): \[\begin{align} p^{(b)}_{mn} &= \frac{b}{d} \sum_{k=0}^{d-1} \sum_{s=0}^{d-1} \omega^{m(k+s)} \frac{1}{1 - \omega^{s+1/2}} \omega^{-kn}= \frac{b}{d} \left( \sum_{k=0}^{d-1} \omega^{k(m-n)} \right) \left( \sum_{s=0}^{d-1} \frac{\omega^{ms}}{1 - \omega^{s+1/2}} \right) \\ &= b \delta_{m,n} \sum_{s=0}^{d-1} \frac{\omega^{ms}}{1 - \omega^{s+1/2}} \;. \end{align}\] The factor \(\delta_{m,n}\) proves that the resulting matrix is diagonal. Let \(x = \omega^{-1/2} = e^{-i\pi/d}\). The diagonal sum becomes: \[T_m = x \sum_{s=0}^{d-1} \frac{(\omega^s)^m}{x - \omega^s} \;.\] Using the root of unity identity \(\sum_{k=0}^{d-1} \frac{z_k^m}{x - z_k} = \frac{d x^{m-1}}{x^d - 1}\) (where \(z_k = \omega^k\)) and noting that \(x^d = (e^{-i\pi/d})^d = -1\), we evaluate \(T_m\) for two cases:

For \(m = 0\): \[T_0 = x \frac{d x^{d-1}}{x^d - 1} = \frac{d x^d}{x^d - 1} = \frac{d(-1)}{-1 - 1} = \frac{d}{2}\] ,

For \(m \neq 0\): \[T_m = x \frac{d x^{m-1}}{x^d - 1} = \frac{d x^m}{-2} = -\frac{d}{2} (e^{-i\pi/d})^m = -\frac{d}{2} \omega^{-m/2} \;.\] Combining \(p^{(a)}_{mn}\) and \(p^{(b)}_{mn}\), the final transformed block is a diagonal matrix with elements: \[p_{mn} = \begin{cases} \alpha & \text{for } m = n = 0 \\ \beta \omega^{-m/2} & \text{for } m = n \neq 0 \\ 0 & \text{for } m \neq n \end{cases}\] where we defined \(\alpha=d(a+b/2)\) (complex) and \(\beta=-d b/2\) (real). In this way, we have identified \(P\), which corresponds to the block \(X^{(12)}\) from 72 expressed in the new basis. The bottom left matrix in 72 is just the conjugate transpose.

Computation of the exponential↩︎

The exact analytical form for the matrix exponential is: \[e^{i X_2} = \begin{pmatrix} \cos(D) & i P D^{-1} \sin(D) \\ iP^\dagger D^{-1} \sin(D) & \cos(D) \end{pmatrix}\] where \(D\) has entries \(d_m = |p_{mm}|\) Because everything inside these blocks is diagonal, we can write down the exact action for any specific \(2 \times 2\) subspace \(m\) (coupling state \(m\) with state \(m+d\)). For a given \(m\), the \(2 \times 2\) submatrix is: \[\exp \left[ i\begin{pmatrix} 0 & p_{mm} \\ p_{mm}^* & 0 \end{pmatrix} \right] = \begin{pmatrix} \cos(|p_{mm}|) & \frac{p_{mm}}{|p_{mm}|} i\sin(|p_{mm}|) \\ i\frac{p_{mm}^*}{|p_{mm}|} \sin(|p_{mm}|) & \cos(|p_{mm}|) \end{pmatrix} \;.\] The subspace corresponding to \(m=0\) is irrelevant because, for \(m=0\), the operator is the identity and the corresponding term in the SOS disappears. Since it is decoupled from the others, we can ignore it. For \(m \neq 0\), we get \[\begin{pmatrix} \cos \beta & i e^{-i \pi m / d} \sin\beta \\ ie^{i \pi m / d} \sin \beta & \cos\beta \end{pmatrix}\] where we assume the interval \(0\le \beta \le \pi/2\). When \(\beta=0\) or \(\beta=\pi/2\) the operators collapse to a combination of operators of the first or second input, respectively. Including a larger interval would just translate into a redefinition of the sign of the operators.

Transformed operators and spectral decomposition↩︎

The final optimal strategy for Bob is identified by \[\begin{align} \widetilde{{B}'}_{j}^{(1)}&=\cos\beta \widetilde{{B}}_j^{(1)}+i \sin \beta \;e^{-i \pi j/d}\widetilde{{B}}_j^{(2)} \\ \widetilde{{B}'}_{j}^{(2)}&=\cos \beta \widetilde{{B}}_j^{(2)}+i \sin \beta \;e^{i \pi j/d}\widetilde{{B}}_j^{(1)} \end{align}\] and we should now find the spectral decomposition of these operators to plug it into the transformed expression of the Bell operator. Recall that \[\widetilde{{B}}_j^{(x)} = \sum_{k=0}^{d-1} \omega^{kj} \ketbra*{\beta_k^{(x)}}{\beta_k^{(x)}} \\ = \frac{1}{d} \sum_{p,q=0}^{d-1} \omega^{(\frac{x}{2}-\frac{1}{4})(p-q)} \ketbra{q}{p} \left( \sum_{k=0}^{d-1} \omega^{k(j-p+q)} \right) \;.\] The sum over \(k\) evaluates to \(d\) if \(p = q + j \pmod{d}\), and \(0\) otherwise. Recalling that \(0 \le j < d\), we have two cases:

• If \(q + j < d\), then \(p = q + j \implies p - q = j\). The phase is \(\omega^{j/4}\).

• If \(q + j \ge d\), then \(p = q + j - d \implies p - q = j - d\). The phase is \(\omega^{(j-d)/4} = -i \omega^{j/4}\).

Let us denote \(q \oplus j \equiv (q + j) \pmod{d}\). Thus, \(\widetilde{{B}}_j^{(1)}\) acts as a cyclic shift with an extra phase: \[\widetilde{{B}}_j^{(1)} \ket{q} = c_q \omega^{j/4} \ket{q \oplus j}\] where \(c_q = 1\) if \(q < d - j\), and \(c_q = -i\) if \(q \ge d - j\). To write a similar expression for \(B_j^{(2)}\), notice that the states \(\ket*{\beta_k^{(1)}}\) and \(\ket*{\beta_k^{(2)}}\) are related by a diagonal phase operator as: \[\ket{\beta_k^{(2)}} = W \ket{\beta_k^{(1)}} \qquad \text{with} \qquad W = \sum_{q} \omega^{q/2} \ketbra{q}{q} = \sum_{q} e^{i \pi q / d} \ketbra{q}{q} \;.\] This implies that \(\widetilde{{B}}_j^{(2)} = W \widetilde{{B}}_j^{(1)} W^\dagger\) and \(\widetilde{{B}}_j^{(2)} \ket{q}\): \[\widetilde{{B}}_j^{(2)} \ket{q} = W \widetilde{{B}}_j^{(1)} e^{-i \pi q / d} \ket{q} = e^{i \pi (q \oplus j)/d} e^{-i \pi q / d} c_q \omega^{j/4} \ket{q \oplus j} = e^{i \pi (q \oplus j - q)/d} \widetilde{{B}}_j^{(1)} \ket{q} \;.\] If \(q + j < d\), the exponent factor is \(e^{i \pi j / d}\), and we get \(e^{i \pi j / d} \widetilde{{B}}_j^{(1)} \ket{q}\). If \(q + j \ge d\), the exponent factor is \(e^{i \pi (j - d) / d} = -e^{i \pi j / d}\), and we get \(-e^{i \pi j / d} \widetilde{{B}}_j^{(1)} \ket{q}\). So, as before, we have a compact form: \[\widetilde{{B}}_j^{(2)} \ket{q} = s_q e^{i \pi j / d} \widetilde{{B}}_j^{(1)} \ket{q}\] where \(s_q = 1\) for \(q < d - j\), and \(s_q = -1\) for \(q \ge d - j\). Substituting this mapping into \(\widetilde{{B}'}_j^{(1)}\), we get: \[\widetilde{{B}'}_j^{(1)} \ket{q} = \left( \cos\beta \widetilde{{B}}_j^{(1)} + i \sin\beta e^{-i \pi j / d} \widetilde{{B}}_j^{(2)} \right) \ket{q} = (\cos\beta + i s_q \sin\beta) \widetilde{{B}}_j^{(1)} \ket{q} = e^{i s_q \beta} \widetilde{{B}}_j^{(1)} \ket{q} \;.\] Similarly, using \(\widetilde{{B}}_j^{(1)} \ket{q} = s_q e^{-i \pi j / d} \widetilde{{B}}_j^{(2)} \ket{q}\), we find: \[\widetilde{{B}'}_j^{(2)} \ket{q} = e^{i s_q \beta} \widetilde{{B}}_j^{(2)} \ket{q} \;.\] Both the final operators are simply the original operators multiplied by a diagonal phase operator. To find the spectral decomposition of \(\widetilde{{B}'}_j^{(1)}\), we set up the eigenvalue equation \(\widetilde{{B}'}_j^{(1)} \ket{w_k} = \lambda_k \ket{w_k}\) and, given its relation to the initial operators, we make the ansatz of a Fourier-like eigenstate \(\ket{w_k} = \frac{1}{\sqrt{d}} \sum_q e^{i \theta_q} \ket{q}\). Applying \(\widetilde{{B}'}_j^{(1)}\) gives: \[\widetilde{{B}'}_j^{(1)} \ket{w_k} = \frac{1}{\sqrt{d}} \sum_{q=0}^{d-1} e^{i \theta_q} e^{i s_q \beta} c_q \omega^{j/4} \ket{q \oplus j} = \lambda_k \frac{1}{\sqrt{d}} \sum_{p=0}^{d-1} e^{i \theta_p} \ket{p} \;.\] To compare the two sides, we need to consider \(p = q \oplus j\) and again we have two regimes. Separating into the two piecewise regimes (\(p \ge j\) and \(p < j\)), we obtain two recurrence conditions for the phases:

• If \(q+j<d\) then \(p=q+j\). Since \(q\ge 0\), it implies that \(p \ge j\). In this case: \[\label{eq:ph1} e^{i \theta_{p-j}} e^{i \beta} \omega^{j/4} = \lambda_k e^{i \theta_p}\tag{73}\]

• If \(q+j \ge d\) then \(p=q+j-d\). Since \(q\le d-1\), it implies \(p < j\). In this case: \[\label{eq:ph2} e^{i \theta_{p-j+d}} e^{-i \beta} (-i) \omega^{j/4} = \lambda_k e^{i \theta_p}\tag{74}\]

To simultaneously solve the two equations 73 and 74 , we assume a linear phase ansatz \(\theta_q = -K q\). Dividing the two equations to isolate \(K\) yields \(e^{-i K d} = i e^{2i \beta}\), which gives the allowed frequencies: \[K = \frac{2\pi}{d} \left( k - \frac{1}{4} - \frac{\beta}{\pi} \right) \quad \text{for } k = 0, 1, \dots, d-1\] Substituting \(K\) back into either recurrence relation yields the exact eigenvalues \(\lambda_k\) and gives the spectral decomposition of \(\widetilde{{B}'}_j^{(1)}\): \[\widetilde{{B}'}_j^{(1)} = \sum_{k=0}^{d-1} \lambda_k \ketbra{u_k(\beta)}{u_k(\beta)}\] where the eigenvectors \(\ket{u_k(\beta)}\) and eigenvalues \(\lambda_k\) are: \[\begin{align} \ket{u_k(\beta)} &= \frac{1}{\sqrt{d}} \sum_{q=0}^{d-1} \exp\left[ -\frac{2\pi i}{d} \left( k - \frac{1}{4} - \frac{\beta}{\pi} \right) q \right] \ket{q} \\ \lambda_k &= \omega^{jk} \exp\left[ i\beta \left( 1 - \frac{2j}{d} \right) \right] \;. \end{align}\] Because \(\widetilde{{B}'}_j^{(2)} = W \widetilde{{B}'}_j^{(1)} W^\dagger\), it shares the exact same eigenvalues \(\lambda_k\), and its eigenvectors are phase-shifted by \(W\): \[\widetilde{{B}'}_j^{(2)} = \sum_{k=0}^{d-1} \lambda_k \ketbra{v_k(\beta)}{v_k(\beta)}\] where the eigenvectors \(\ket{v_k(\beta)} = W \ket{u_k(\beta)}\) are given by: \[\begin{align} \ket{v_k(\beta)} = \frac{1}{\sqrt{d}} \sum_{q=0}^{d-1} \exp\left[ -\frac{2\pi i}{d} \left( k - \frac{3}{4} - \frac{\beta}{\pi} \right) q \right] \ket{q} \;. \end{align}\] As a quick check, setting \(\beta = 0\) perfectly collapses \(\ket*{u_k(\beta)}\) into \(\ket*{\beta_k^{(1)}}\), \(\ket*{v_k(\beta)}\) into \(\ket*{\beta_k^{(2)}}\), and \(\lambda_k\) into \(\omega^{jk}\), recovering the original operators \(\widetilde{{B}}_j^{(1)}\) and \(\widetilde{{B}}_j^{(2)}\). The measurements for the standard SATWAP inequality are recovered by choosing \(\beta=\pi/4\).

References↩︎

---

1. lorenzo.coccia@unipd.it↩︎

2. Indeed, because the transformed operators are built from the initial elements, the space spanned by the new operators is contained within the space spanned by the original operators. At the same time, since the transformation is invertible, the space spanned by the original operators is contained entirely within the space spanned by the new operators. Therefore, the two spaces must be identical.↩︎