\section{Correlation form of steering functionals}
\label{AppA}

Let us consider the standard form of the steering functional, i.e., 
\begin{equation}\label{SteeringFunctional}
    \sum_{b=0}^{d-1}\sum_{y=1}^N \Tr\left(F_{b|y}\sigma_{b|y}\right),
\end{equation}
where $\sigma_{b|y}$ form the assemblage and $F_{b|y}$ are some positive
semi-definite operators acting on $\mathbbm{C}^d$. Let us then consider a particular class of such functionals in which for every $y=1,\ldots,N$, the operators $F_{a|x}$ sum up to the identity, that is, 
\begin{equation}
    \sum_{a}F_{a|x}=\mathbbm{1},
\end{equation}
and therefore form a quantum measurement. With the aid of Fourier transform
we can represent each such a measurement by another set of operators 
\begin{equation}\label{Akx}
    A_{k|x}=\sum_{a=0}^{d-1}\omega^{ak} F_{a|x}
\end{equation}
with $k=0,\ldots,d-1$ and $\omega=\mathrm{exp}(2\pi\mathbbm{i}/d)$, which for a moment we can call generalized observables; in fact if $F_{a|x}$ forms a projective measurement, then $A_{k|x}$ are unitary observables with outcomes labelled by the powers of $\omega$ (see also below). 

The inverse transform gives
\begin{equation}
F_{a|x}=\frac{1}{d}\sum_{k=0}^{d-1}\omega^{-ak} A_{k|x}.   
\end{equation}




Clearly, we can use the same approach to represent Bob's measurements $N_y=\{N_{b|y}\}_b$ in terms of generalized observables, that is, 
\begin{equation}
    B_{k|y}=\sum_{b=0}^{d-1}\omega^{bk} N_{b|y},
\end{equation}
and
\begin{equation}
    N_{b|y}=\frac{1}{d}\sum_{k=0}^{d-1}\omega^{-kb}B_{k|y}.
\end{equation}
It is worth recalling here that while the dimension of Alice's observables is known to be $d$, Bob's measurements act on a Hilbert space of in principle unknown but finite dimension. 

Let us now go back to the steering functionals (\ref{SteeringFunctional}). By exploiting the fact that $\sigma_{b|y}=\Tr_B[(\mathbbm{1}_A\otimes N_{b|y})\rho_{AB}]$, it can be rewritten as a combination of 
expectation values of the above generalized observables, 
\begin{equation}\label{SteeringFunctional2}
    \sum_{b=0}^{d-1}\sum_{y=1}^N \Tr\left(F_{b|y}\sigma_{b|y}\right)=\frac{1}{d}
   \sum_{k=0}^{d-1} \sum_{y=1}^N\left\langle A_{k|y}\otimes B_{-k|y}\right\rangle.
\end{equation}

We then further restrict the class of the steering functionals to the case when the operators $F_{a|x}$ and $N_{b|y}$ form projective measurements for all $x$ and $y$, meaning that they satisfy 
\begin{equation}
    F_{b|y}F_{b'|y}=\delta_{b,b'}F_{b|y},\qquad
    N_{b|y}N_{b'|y}=\delta_{b,b'}N_{b|y}.
\end{equation}
It is not difficult to see that under the above assumption, $A_{k|y}$ as well as $B_{k|y}$ are all unitary operators with eigenvalues $\omega^{bk}$ $(b=0,\ldots,d-1)$, and, moreover, $A_{k|y}$ and $B_{k|y}$ for any $y$ are simply the $k$th powers of $A_{1|y}$ and $B_{1|y}$, respectively. Using this fact,  denoting
\begin{equation}
    A_{1|y}\equiv A_{y}, \qquad B_{1|y}=B_{y},
\end{equation}
and dropping the $1/d$, 
we finally arrive at our steering functional of the form
\begin{equation}\label{SteeringFunctional21}
    \langle \mathcal{B}_d\rangle=
   \sum_{k=1}^{d-1} \sum_{y=1}^N\left\langle A_{y}^k\otimes B_{y}^{d-k}\right\rangle,
\end{equation}
where we also got rid of the term for $k=0$ because 
the expectation value for this term is simply one.
Let us also notice that without any loss of generality we
can replace in Eq. (\ref{SteeringFunctional21}) $d-k$ by $k$
as such an operation corresponds to a simple permutation of the 
eigenvalues of $B_i$, which leads us to functionals considered in the 
main text. 

Let us stress here that this the above form of the steering functional is particularly useful as far as self-testing is concerned because
it is direct to find its maximal quantum value, i.e., $\beta_Q=N(d-1)$ as well as a quantum realization achieving it. Moreover, for any quantum realization achieving $\beta_Q$, the following relations 
\begin{equation}
    A_y^k\otimes B_y^{-k}\ket{\psi}=\ket{\psi}
\end{equation}
must hold for any $y$ and $k$.




\section{Bounding the maximal classical value of our functionals}
\label{AppB}


In this appendix we first derive an upper bound on the maximal classical value
of the steering expression \eqnref{classBoun}, and show that for any set of 
genuinely $d$-outcome incompatible measurements the corresponding steering inequality is nontrivial, that is, $\beta_L<\beta_Q$. We then provide a self-testing proof in the exact case and in the case when, due to noises or experimental imperfections, the violation of the steering functional is suboptimal. 

For further convenience let us recall here that our 
steering functionals read
\begin{eqnarray}\label{SteeOp}
\langle\mathcal{B}_d\rangle=\sum_{i=1}^{N}\sum_{k=1}^{d-1} \left\langle A_i^{k}\otimes B_i^{k}\right\rangle,
\end{eqnarray}
where, as compared to Eq. (\ref{SteeringFunctional21}) we have 
replaced $B_i^{d-k}$ by $B_i^k$ which corresponds to a simple permutation
of the outcomes of $B_i$. Their maximal quantum values are $\beta_{Q}=N(d-1)$. 

\subsection{Classical bound}

Let us now demonstrate how to derive the bound (\ref{classBoun}). To this end, we consider assume that the assemblage $\{\sigma_{b|y}\}$ admit a local hidden state model, that is, 
\begin{equation}
    \sigma_{b|y}=\sum_{\lambda}p(\lambda)p(b|y,\lambda)\rho_{\lambda},
\end{equation}
where $\lambda$ is some random variable with some probability distribution $p(\lambda)$, 
$p(b|y,\lambda)$ a probability distributions according to which Bob returns the outcomes $a$ given the measurement choice $y$ and the hidden variable $\lambda$, and $\rho_{\lambda}$
are some states acting on Alice's Hilbert space. In terms of the joint probability 
distributions such an LHS model expresses as
\begin{eqnarray}
p(a,b|x,y)=\sum_{\lambda}p(\lambda)p(a|x,\rho_\lambda)p(b|y,\lambda),
\end{eqnarray}
where $p(a|x,\rho_\lambda)$ is the probability of obtaining the outcome $a$
when measuring the observable $A_x$ on the state $\rho_{\lambda}$. Then, in terms of the expectation values appearing in (\ref{SteeringFunctional2}), 
the above rewrites as 
\begin{equation}
\left\langle A_{k|x}\otimes B_{l|y}\right\rangle_{\mathrm{LHS}}=
\sum_{\lambda}p(\lambda)\langle A_{k|x}\rangle_{\rho_{\lambda}} 
\langle B_{l|y}\rangle_{\lambda},
\end{equation}
where $\langle A_{k|x}\rangle_{\rho_{\lambda}}=\Tr[A_{k|x}\rho_{\lambda}]$ with $A_{k|x}$ being operators representing the $x$th measurement performed by Alice 
[cf. Eq. (\ref{Akx})]. Then, $\langle B_{l|y}\rangle_{\lambda}$ stands for the Fourier transform of $p(b|y,\lambda)$, i.e., 
\begin{equation}
   \langle B_{l|y}\rangle_{\lambda}=\sum_{b}\omega^{bl}p(b|y,\lambda).
\end{equation}

Now, our steering expression (\ref{SteeringFunctional21}) for LHS models reads
\begin{equation}
\langle\mathcal{\hat{B}}_d\rangle_{\mathrm{LHS}}=\sum_{\lambda}p(\lambda) \sum_{y=1}^N\sum_{k=1}^{d-1}\langle A_{k|y}\rangle_{\rho_{\lambda}} 
\langle B_{-k|y}\rangle_{\lambda},
\end{equation}
and by noting that $|\langle B_{l|y}\rangle_{\lambda}|\leq 1$ for any choice of $y$, $l$ and $\lambda$, its absolute value can be upper bounded as 
\begin{equation}
|\langle\mathcal{\hat{B}}_d\rangle_{\mathrm{LHS}}|\leq \sum_{\lambda}p(\lambda) \sum_{i=1}^N\sum_{k=1}^{d-1}|\langle A_i^k\rangle_{\rho_{\lambda}}|,
\end{equation}
and further as 
\begin{equation}\label{dupa}
|\langle\mathcal{\hat{B}}_d\rangle_{\mathrm{LHS}}|\leq \max_{\rho} \sum_{i=1}^N\sum_{k=1}^{d-1}|\langle A_i^k\rangle_{\rho}|,
\end{equation}
where we took advantage of the fact that the right-hand side of 
the above inequality is a convex combination and that by optimizing it over states $\rho$ acting on Alice's Hilbert space we cannot decrease its value. We thus obtain the bound (\ref{classBoun}).

A direct, but important for our considerations corollary that stems from the above bound is that for any set of genuinely incompatible measurements 
$A_i$, the right hand side of (\ref{dupa}) is smaller than $N(d-1)$, and thus $\beta_L<N(d-1)$, meaning that for any such a set our steering inequality is nontrivial.  
To see this explicitly let us assume that 
\begin{equation}
    \sum_{i=1}^N\sum_{k=0}^{d-1}|\langle A_i^k\rangle_{\rho}|=N(d-1),
\end{equation}
which implies that $|\langle A_i^k\rangle_{\rho}|=1$ for any $i$, $k$ and some $\ket{\psi}$. Due to the fact that $A_i^k$ are unitary operators, this means that $\ket{\psi}$ is their eigenvector, possibly with different eigenvalues. Consequently, these observables share 
a common subspace spanned by this vectors which contradicts that they are genuinely incompatible.



\section{Self-testing}

Here we present proofs of Lemmas \ref{lemma4} and \ref{Lem2}, and Theorems \ref{Theo1M} and \ref{Thm222}.

\label{AppC}

\setcounter{lem}{0}

\begin{lem}
Consider a set of $N$ genuinely $d$-outcome incompatible observables $A_i$ on $\mathbbm{C}^d$. If for some nonzero $P\geq 0$ acting on $\mathbbm{C}^d$, $[P,A_i]=0$ for every $i=1,\ldots,N$, then $P=\mu \mathbbm{1}$ with $\mathbbm{1}$ being the identity acting on $\mathbbm{C}^d$ and $\mu> 0$.
\end{lem}
\begin{proof}
The proof is by contradiction. Let us assume that the matrix $P$ is not proportional to the identity. Since it is positive semi-definite, 
we can decompose it as
\begin{equation}
    P=\sum_{i}\lambda_i P_i,
\end{equation}
where $\lambda_i\geq 0$ are the eigenvalues of $P$ which without any loss of generality we can assume to be all different, and $P_i$ are the corresponding orthogonal eigenprojections. 

We first show that $\mathrm{rank}(P)=d$, that is, that $\lambda_i>0$ for all $i$. To this end, let us assume that $P$ has a nonempty kernel, that is there exists
$\ket{\psi}\in\mathbbm{C}^d$ such that $P\ket{\psi}=0$. Then, we can always write
\begin{equation}
    A_i\ket{\psi}=\alpha_i\ket{\psi}+\beta_i\ket{\psi_i},
\end{equation}
where $\alpha_i,\beta_i\in\mathbbm{C}^d$ and $\ket{\psi_i}$ are some 
vectors such that $\langle \phi_i|\psi\rangle=0$ for any $i$. Applying now $P$ 
to both sides of the above, we obtain $PA_i\ket{\psi}=\beta_i\ket{\psi_i}$. 
Then, due to the fact that $P$ commutes with all $A_i$, we
have $\beta_i\ket{\psi_i}=PA_i\ket{\psi}=A_iP\ket{\psi}=0$, that is, 
$\beta_i\ket{\psi_i}$ for all $i$. This finally implies that $\ket{\psi}$ is an 
eigenvector of all observables $A_i$ which contradicts the assumption that 
these observables are genuinely incompatible.

Now, it is a known fact that if $[P,A_i]=0$, then $A_i$ is of block form 
\begin{equation}\label{dupa12}
 A_i=A_i^{(1)}\oplus \ldots \oplus A_i^{(m)},
\end{equation}
where each block $A_i^{(j)}$ acts on $\mathrm{supp}(P_j)$ and $m$ is the number of different eigenvalues of $m$. For completeness, let us add here a short proof of this fact. To this aim, we rewrite $[P,A_i]=0$ as
\begin{equation}
    A_iP=PA_i.
\end{equation}
We then sandwich this equation with $P_m$ and $P_n$ for $m\neq n$, which 
gives
\begin{equation}
    P_mA_iPP_n=P_mPA_iP_n.
\end{equation}
Noting that $PP_n=\lambda_nP_n$ and $P_mP=\lambda_m P_m$, 
we can rewrite the above as
\begin{equation}
    (\lambda_n-\lambda_m)P_mA_i P_n=0,
\end{equation}
and since $\lambda_n\neq \lambda_m$ for $n\neq m$, we conclude that 
$P_mA_i P_n=0$ for any pair $m\neq n$. Of course, this conclusion can be drawn
for each $A_i$, and thus we obtain (\ref{dupa12}). 

The fact that all observables can be written as in (\ref{dupa12}) implies that 
$\mathrm{supp}(P_i)$ is an invariant subspace of every $A_i$ of dimension lower than $d$, which contradicts the assumption. This completes the proof.
\end{proof} 

Let us now provide a full proof of our main result.

\setcounter{thm}{0}

\begin{thm}\label{Theo1} 
Assume that the steering inequality \eqref{SteeIn} is maximally violated by a state $\ket{\psi}_{AB}\in\mathbbm{C}^d\otimes\mathcal{H}_B$ and observables $A_i,B_i \ (i\in \{1,2,\ldots,N\})$ acting on, respectively, $\mathbbm{C}^d$ and $\mathcal{H}_B$ such that the observables $A_i$ on the trusted side are $d$-outcome genuinely incompatible. Then, the following statement holds true for any $d$: there exist a local unitary transformation on Bob side, $U_B:\mathcal{H}_B\rightarrow\mathcal{H}_B$, such that
\begin{eqnarray}\label{Theo1.1}
\forall i, \quad U_B\overline{B}_iU_B^{\dagger}=A_i^{*}
\end{eqnarray}
where $\overline{B}_i$ is $B_i$ projected onto the support of Bob's state $\rho_B\in\mathbbm{C}_d$ and 
\begin{eqnarray}\label{Theo1.2}
(\mathbbm{1}_A\otimes U_B)\ket{\psi}_{AB}=
\ket{\phi_d^+}.
\end{eqnarray}
\end{thm}
\begin{proof}
The optimal quantum value $\beta_{Q}=N(d-1)$ of the steering inequality \eqref{SteeIn} is the same as its algebraic value. Thus, we have the following set of relations
\begin{eqnarray}\label{SOS}
\quad  A_i^{k}\otimes B_i^{k} \ket{\psi}_{AB} = \ket{\psi}_{AB}.
\end{eqnarray}
for all $i$ and $k$. For clarity the rest of the proof is divided into 
three main steps.\\

\textit{Step 1.} Let us first use these relations
to prove that the fact that $A_i$ form a set of
genuinely incompatible measurements implies that $\mathrm{rank}(\rho_A)=d$
with $\rho_A$ being the reduced density matrix corresponding to Alice;
in other words, any state maximally violating our inequality 
for genuinely incompatible measurements must of maximal local 
rank. To this end, we assume that $\mathrm{rank}(\rho_A)<d$ and consider (\ref{SOS}) 
for $k=1$, which implies that  
\begin{eqnarray}\label{SOS22}
\quad  \Pi_A\, A_i\,\Pi_A\otimes B_i \ket{\psi}_{AB} = \ket{\psi}_{AB},
\end{eqnarray}
where $\Pi_A$ is the projector onto the support of 
$\rho_A$. Denoting then $\overline{A}_i\equiv \Pi_A\, A_i\,\Pi_A$
and owing to the fact that $B_i$ are unitary and satisfy 
$B_i^d=\mathbbm{1}$, we conclude that $\overline{A}_i$ are unitary 
matrices on $\mathrm{supp}(\rho_A)$ and satisfy $\overline{A}_i^d=\Pi_A$.
It is then not difficult to see that since $A_i$ are unitary, they 
must be of block form 
\begin{equation}
    A_i=\overline{A}_i\oplus A'_i 
\end{equation}
for some square unitary matrix $A'_i$ of dimension $[d-\mathrm{rank}(\rho_A)]\times[d-\mathrm{rank}(\rho_A)]$.
However, this implies that all the observables $A_i$ have a
common invariant subspace of dimension lower $d$, contradicting the
assumption that they are genuinely incompatible on $\mathbbm{C}^d$. 
As a consequence, $\ket{\psi}_{AB}$ is of full local rank, $\mathrm{rank}(\rho_A)=d$ as well as $\mathrm{rank}(\rho_B)=d$.\\

\textit{Step 2.} Using a similar reasoning we now show that Bob's observables 
all split into a direct sum 
\begin{equation}\label{block}
    B_i=\overline{B}_i\oplus B_{i}^{\mathrm{rest}},
\end{equation}
where $\overline{B}_i=\Pi_B\,B_i\,\Pi_B$ are $d\times d$ matrices acting 
on the support of $\rho_B$, whereas $B_i^{\mathrm{rest}}$ are of unknown dimension and are defined outside $\mathrm{supp}(\rho_B)$; in particular, $B_i^{\mathrm{rest}}\ket{\psi}_{AB}=0$. Both $\overline{B}_i$ and $B_i^{\mathrm{rest}}$ are unitary and their spectra belong to $\{1,\ldots,\omega^{d-1}\}$. 

To prove the above statement, let us first notice that 
each $B_i$ can be written in a block form,
\begin{equation}
    B_i=\left(
    \begin{array}{cc}
    \overline{B}_i & C_i\\
    D_i & B_i^{\mathrm{rest}}
    \end{array}
    \right),
\end{equation}
where, as before, $\overline{B}_i=\Pi_B\,B_i\,\Pi_B$
and $\overline{B}_i^{\mathrm{rest}}=\Pi_B^{\perp}\,B_i\,\Pi_B^{\perp}$
where $\Pi_B$ and $\Pi_B^{\perp}$ are the projectors onto the support of $\rho_B$ and its in $\mathcal{H}_B$, respectively. Let us then go back to Eq. (\ref{SOS}) for $k=1$ and project it onto the support of $\rho_B$, that is, 
\begin{equation}\label{Jumila}
 A_i\otimes \Pi_B\, B_i\,\Pi_B \ket{\psi}_{AB} \equiv A_i\otimes  \overline{B}_i\ket{\psi}_{AB} = \ket{\psi}_{AB}.
\end{equation}
By applying $A_i^{\dagger}\otimes \overline{B}_i^{\dagger}$ to this identity and using the fact that $A_i$ are unitary, we conclude that $\overline{B}_i^{\dagger}\overline{B}_i=\mathbbm{1}_d$. In a fully analogous way 
we obtain from Eq. (\ref{SOS}) for $k=d-1$ that $\overline{B}_i\overline{B}_i^{\dagger}=\mathbbm{1}_d$ and thus 
all $\overline{B}_i$ are unitary. On the other hand, by applying $A_i\otimes\overline{B}_i$ to Eq. (\ref{Jumila}) $d-1$ times, we
obtain that $\overline{B}_i^d=\mathbbm{1}_d$ and thus their spectra
belong to $\{1,\ldots,\omega^{d-1}\}$. 

It is then direct to see that since both $B_i$ and $\overline{B}_i$ 
are unitary on the respective Hilbert spaces, the off-diagonal blocks in Eq. (\ref{Jumila}) must vanish, $C_i=D_i=0$, leading to the desired block form 
(\ref{block}). This finally implies that $B_i^{\mathrm{rest}}$ must also be
unitary and their eigenvalues are powers $\omega$. \\




\textit{Step 3.} Let us then consider the Schmidt decomposition of 
$\ket{\psi_{AB}}$,
\begin{equation}
    \ket{\psi}_{AB}=\sum_{i=0}^{d-1}\lambda_i\ket{e_i}\ket{f_i},
\end{equation}
where due to the fact that $\mathrm{rank}(\rho_A)=d$, the Schmidt coefficient $\lambda_i$ $(i=0,\ldots,d-1)$ are all positive and satisfy $\sum_{i}\lambda_i^2=1$. Moreover, the local vectors $\ket{e_i}$ and $\ket{f_i}$ 
form orthonormal bases in $\mathbbm{C}^d$ and $\mathcal{H}_B$, respectively.


Now, combining Eqs. \eqref{PW} and \eqref{SOS00} for $k=1$, we obtain
\begin{eqnarray}
(A_iP_A\otimes \widetilde{B}_i)\ket{\phi^{+}_d}=(P_A\otimes \mathbbm{1}_B)\ket{\phi^{+}_d}\nonumber,
\end{eqnarray}
where $\widetilde{B}_i=U_B \overline{B}_i U_B^\dagger$.
Using then the fact that $R\otimes Q\ket{\phi_d^{+}}=RQ^T\otimes\mathbbm{1}\ket{\phi_d^{+}}$
holds for any two $d\times d$ matrices $R$ and $Q$,
and omitting the state $\ket{\phi_d^+}$, we obtain 
a set of matrix equation 
\begin{equation}\label{CULO-A}
A_iP_A\widetilde{B}_i^T=P_A.
\end{equation}
After multiplying it by its hermitian conjugation, and 
rearranging some terms, we finally arrive at a simple condition 
for the matrix $P_A$, $[A_i,P_A^2]=0$, which due to the fact that 
$P_A\geq0$ is equivalent to
\begin{equation}
    [A_i,P_A]=0
\end{equation}
for all $i$. Due to the facts that $P_A\geq 0$ and that $A_i$
are genuinely incompatible, the above lemma tell us that $P_A$
is proportional to the identity, proving Eq. (\ref{Theo1.2M}).
Then, plugging this particular form of $P_A$ into Eq. (\ref{CULO-A})
we obtain (\ref{Theo1.1M}), which completes the proof.
\end{proof}




Now, we move on to prove that any set of observables whose eigenvectors
form mutually unbiased bases are genuinely incompatible. 
\begin{lem}\label{Lem2A}
Any two $d$-outcome observables whose eigenbases are mutually unbiased are genuinely incompatible.
\end{lem}
\begin{proof}
The proof is by contradiction. Let us begin by taking two $d-$dimensional mutually unbiased bases denoted by, $\{\ket{s_i}\}$ and $\{\ket{t_j}\}$ such that $|\langle s_i|t_j\rangle|^2=1/d$ for all $i,j\in\{0,1,\ldots,d-1\}$. Let us then construct observables from these bases, 
\begin{equation}
A_1=\sum_{i=0}^{d-1}\omega^{f(i)}\ket{s_i}\!\bra{s_i},\qquad  A_2=\sum_{i=0}^{d-1}\omega^{g(i)}\ket{t_i}\!\bra{t_i},
\end{equation}
where $f$ and $g$ are some permutations of the $d$-element set $\{0,\ldots,d-1\}$,
and assume that they share a common invariant subspace $V\subsetneq \mathbbm{C}^d$ such that $\dim V<d$. This means that both can be written as 
$A_i=A_i'\oplus A_i''$, where $A_i'$ acts on $V$, whereas 
$A_i''$ on its complement in $\mathbbm{C}^d$. Then, $A_1$ has eigenvectors
that are orthogonal to eigenvectors of $A_2$, which contradicts that 
the eigenvectors of $A_1$ and $A_2$ are mutually unbiased.
\end{proof}
For a note, we construct a set of genuinely incompatible measurements for $d=4$ which are not pairwise genuinely incompatible. For this, we would require the underlying lemma.
\begin{lem}\label{lemma3}
Two invertible and diagonalizable matrices $T$ and $T'$ of dimension $D$, share a common non-trivial invariant subspace of dimension $d<D$ iff there exist $d$ eigenvectors of $T'$ which can be written as a linear combination of $d$ eigenvectors of $T$.
\end{lem}  

\begin{proof}
Let us first recall that any non-trivial invariant subspace of a matrix $T$ is spanned by its eigenvectors which implies that if $V_d$ is a $d-$dimensional common invariant subspace of $T$ and $T'$, then $V_d$ is spanned by $d$ eigenvectors of $T$ as well as $d$ eigenvectors of $T'$. Since, eigenvectors of any matrix are linearly independent, this implies that there exist $d$ eigenvectors of $T'$ which can be written as linear combinations of $d$ eigenvectors of $T$. 

Also, if $d$ eigenvectors of $T'$, denoted by $\{\ket{e'_i}\}|_{i=1,\ldots,d}$ can be written as a linear combination of $d$ eigenvectors of $T$, denoted by $\{\ket{e_i}\}|_{i=1,\ldots,d}$, then there exist a common invariant subspace $V_d$ of $T$ and $T'$ which is spanned by $\{\ket{e'_i}\}|_{i=1,\ldots,d}$. This completes the proof.
\end{proof}  

The example is stated below.
\begin{small}
\begin{eqnarray}\label{C19}
    A_1=\frac{1}{2}\begin{pmatrix}
1+\mathbbm{i} & 1-\mathbbm{i} & 0 & 0\\
1-\mathbbm{i} & 1+\mathbbm{i} & 0 & 0\\
0 & 0 & -2 & 0\\
0 & 0 & 0 & -2\mathbbm{i}
\end{pmatrix}, \quad   
A_2=\frac{1}{2}\begin{pmatrix}
2 & 0 & 0 & 0\\0 & -1+\mathbbm{i} & 1+\mathbbm{i} & 0\\0 &
1+\mathbbm{i}& -1+\mathbbm{i} & 0\\
0 & 0 & 0 & -2\mathbbm{i}
\end{pmatrix}, \quad
 A_3=\frac{1}{2}\begin{pmatrix}
2 & 0 & 0 & 0\\0 & 2\mathbbm{i} & 0 & 0\\0 & 0 &-1-\mathbbm{i} & \mathbbm{i}-1 \\0 &
0 & \mathbbm{i}-1& -1-\mathbbm{i} \end{pmatrix}
\end{eqnarray}
\end{small}

As discussed before, any non-trivial invariant subspace of a matrix $T$ is spanned by its eigenvectors of $T$. Using this fact, any non-trivial invariant subspace of $A_1$ is spanned by the vectors $\{\ket{+_{01}},\ket{-_{01}},\ket{2},\ket{3}\}$, any non-trivial invariant subspace of $A_2$ is spanned by the vectors $\{\ket{0},\ket{+_{12}},\ket{-_{12}},\ket{3}\}$ 
and any non-trivial invariant subspace of $A_3$ is spanned by the vectors $\{\ket{0},\ket{1},\ket{+_{23}},\ket{-_{23}}\}$ where $|\pm_{ij}\rangle = (|i\rangle \pm |j\rangle)/\sqrt{2}$.

Using Lemma \ref{lemma3}, we can conclude that $A_1$ and $A_3$ have two common nontrivial invariant subspaces, spanned by $\{\ket{0},\ket{1}\}$ and $\{\ket{2},\ket{3}\}$. However, these two are not invariant subspaces of $A_2$ because it is not possible to represent an arbitrary vector from those subspaces as linear combinations of two of the eigenvectors of $A_2$. Thus, we can conclude that there is no common invariant subspace shared between $A_1, A_2$ and $A_3$.    




\subsection{Robust self-testing}

\begin{thm}
Consider a pair of genuinely $d-$outcome incompatible observables $A_1=X_dZ_d^{l}$ with $l=0,\ldots,d-1$ and $A_2=Z_d$, and assume that the following steering inequality 
\begin{eqnarray}\label{InexIn}
\sum_{i=1}^{2}\sum_{k=1}^{d-1}\left\langle A_i^{k}\otimes B_i^{k}\right\rangle_{\ket{\psi}}\geq 2(d-1) -\varepsilon
\end{eqnarray}
is violated $\epsilon$-close to the maximal violation by a state 
$\ket{\psi}\in\mathcal{H}_{A}\otimes\mathcal{H}_B$ and observables $B_i \ (i=1,2)$ acting on $\mathcal{H}_B$. Then, for any $d\geq 2$ there exists a unitary operation $U_B:\mathcal{H}_B\to \mathcal{H}_B$ such that
\begin{eqnarray}\label{Rob1}
\left\|(\mathbbm{1}\otimes U_B)(\mathbbm{1}\otimes B_i^{k})\ket{\psi}-\mathbbm{1}\otimes (B'_i)^{k}\ket{\phi^+_d}\right\|\leq \sqrt{2\varepsilon}+2\sqrt{d}\sqrt[4]{2\varepsilon}
\end{eqnarray}
and
\begin{eqnarray}\label{Rob2}
\left\|B_i^{k}-(B'_i)^{k}\right\|_2^2\leq \sqrt{d}\sqrt{2\varepsilon}
\left(1+4\sqrt{d}\sqrt{2\varepsilon}\right)
\end{eqnarray}
with $k=0,\ldots,d-1$, where $B'_i$ are the ideal observables of Bob, i.e., $B_i'=A_i^*$, and $\|\cdot\|_2$ stands for the Hilbert-Schmidt norm.
\end{thm}
\begin{proof}The first step in proving the above statements 
is to observe that 
(\ref{InexIn}) implies the following set of inequalities
\begin{eqnarray}\label{App:1}
\left|\langle\psi| A_i^k\otimes B_i^k|\psi\rangle\right|\geq \mathrm{Re}\left(\langle\psi| A_i^k\otimes B_i^k|\psi\rangle\right)\geq 1-\varepsilon
\end{eqnarray}
for any $i=1,2$ and $k=1,\ldots,d-1$. 

The second observation we will exploit in our proof is that the shared state $\ket{\psi}$ can always be written in the basis of $Z_d$ as
\begin{eqnarray}\label{stateRob}
\ket{\psi}=\sum_{i=0}^{d-1}\alpha_i\ket{i}\ket{b_i},
\end{eqnarray}
where $\ket{b_i}$ are some, not necessarily orthogonal vectors from $\mathcal{H}_B$, and $\alpha_i$ are non-negative numbers satisfying $\sum_{i=0}^{d-1}\alpha_i^2=1$.

Let us first exploit Ineq. (\ref{App:1}) for $i=1$. Exploiting the fact that 
$(X_dZ_d^l)^{k}\ket{i}=\omega^{kl\left(i+\frac{k-1}{2}\right)}\ket{i+k}$, where the addition is mod $d$, we obtain
\begin{equation}
\sum_{i}\alpha_i\alpha_{i+k}\mathrm{Re}\left(\omega^{kl\left(i+\frac{k-1}{2}\right)}\langle b_{i+k}|B_1^k|b_i\rangle\right)\geq 1-\varepsilon
\end{equation}
which by virtue of the facts that $\mathrm{Re}(z)\leq |z|$
for any $z\in\mathbbm{C}$ and that $\left|\omega^{kl\left(i+\frac{k-1}{2}\right)}\langle b_{i+k}|B_1^k|b_i\rangle\right|\leq 1$ because $B_1$ is unitary and $\ket{b_i}$ are normalized, implies
\begin{equation}
\sum_{i}\alpha_i\alpha_{i+k}\geq 1-\varepsilon
\end{equation}
for any $k=0,\ldots,d-1$. By summing the above relation over $k$, one obtains
\begin{equation}\label{App:3}
\sum_{i,k=0}^{d-1}\alpha_i\alpha_{i+k}=\left(\sum_{i=0}^{d-1}\alpha_{i}\right)^2\geq d(1-\varepsilon)
\end{equation}
which gives 
\begin{equation}\label{App:2}
\sum_{i=0}^{d-1}\alpha_{i}\geq \sqrt{d}\sqrt{1-\varepsilon}.
\end{equation}
We can use this inequality to show that each $\alpha_i$ is 
"$\varepsilon$-close" to 
$1/\sqrt{d}$. To this end, let us write that
\begin{eqnarray}
    \sum_{i=0}^{d-1}\left(\alpha_i-\frac{1}{\sqrt{d}}\right)^2&=&\sum_{i=0}^{d-1}
    \alpha_i^2-\frac{2}{\sqrt{d}}\sum_{i=0}^{d-1}\alpha_i+1\nonumber\\
    &=&2\left(1-\frac{1}{\sqrt{d}}\sum_{i=0}^{d-1}\alpha_i\right),
\end{eqnarray}
where we used that $\sum_{i}\alpha_i^2=1$. Now, by using (\ref{App:2})
we can upper bound the above expression as
\begin{eqnarray}
    \sum_{i=0}^{d-1}\left(\alpha_i-\frac{1}{\sqrt{d}}\right)^2&\leq &
    2(1-\sqrt{1-\varepsilon})\nonumber\\
    &\leq& 2\varepsilon,
\end{eqnarray}
where we used $1-\sqrt{1-\varepsilon}\leq \varepsilon$, which in turn is a simple consequence of the fact that $\sqrt{1-\varepsilon}\geq 1-\varepsilon$ for any $0\leq\varepsilon\leq1$. This finally implies that for any $i$,
\begin{equation}\label{App:4.1}
\frac{1}{\sqrt{d}}-\sqrt{2\varepsilon}\leq \alpha_i\leq 
\frac{1}{\sqrt{d}}+\sqrt{2\varepsilon}.
\end{equation}

In a similar manner we can prove that 
\begin{equation}\label{App:4}
\frac{1}{d}-\sqrt{2\varepsilon}\leq \alpha_i\alpha_{i+j}\leq 
\frac{1}{d}+\sqrt{2\varepsilon}.
    \end{equation}
for any $i,j=0,\ldots,d-1$. Precisely, let us consider the following expression
\begin{eqnarray}
    \sum_{i,j=0}^{d-1}\left(\alpha_{i}\alpha_{i+j}-\frac{1}{d}\right)^2&=&2-\frac{2}{d}\sum_{i,j=0}^{d-1}\alpha_i\alpha_{i+j}\nonumber\\
    &\leq &2-\frac{2}{d}d(1-\varepsilon)=2\varepsilon,
\end{eqnarray}
where we have used Eq. (\ref{App:3}). By noting that each term in the above sum must also be upper bounded by $2\varepsilon$ we arrive at (\ref{App:4}).

Let us now use conditions (\ref{App:1}) for $A_2=Z_{d}$, which by using the fact that $Z_d^k\ket{i}=\omega^{ik}\ket{i}$ can be rewritten 
as
\begin{equation}
    \sum_{i=0}^{d-1}\alpha_i^2\,\mathrm{Re}\left(\omega^{ik}\langle b_i|B_2^k|b_i\rangle\right)\geq 1-\varepsilon.
\end{equation}
Using then the fact that $\mathrm{Re}(\omega^{ik}\langle b_i|B_2^k|b_i\rangle)\leq 1$ for any $i,k$, we have
\begin{equation}
    \sum_{i\neq j}\alpha_i^2+\alpha_j^2\,\mathrm{Re}\left(\omega^{jk}\langle b_j|B_2^k|b_j\rangle\right)\geq 1-\varepsilon,
\end{equation}
which by virtue of $\sum_i\alpha_i^2=1$ implies
\begin{equation}
    \alpha_j^2\left[1-\mathrm{Re}\left(\omega^{jk}\langle b_j|B_2^k|b_j\rangle\right)\right]\leq \varepsilon.
\end{equation}
Using then the left-hand side inequality of Eq. (\ref{App:4.1}) 
and the fact the term in the square brackets is non-negative, we obtain
\begin{equation}
    \frac{1}{d}\left[1-\mathrm{Re}\left(\omega^{jk}\langle b_j|B_2^k|b_j\rangle\right)\right]\leq \varepsilon+\sqrt{2\varepsilon}\leq 2\sqrt{2\varepsilon},
\end{equation}
which, by rearranging terms gives
\begin{equation}\label{App:5}
\mathrm{Re}\left(\omega^{jk}\langle b_j|B_2^k|b_j\rangle\right)\geq 1-2d\sqrt{2\varepsilon}.
\end{equation}
Let us now consider the eigendecompositon of $B_2$,
\begin{equation}
    B_2=\sum_{i=0}^{d-1}\omega^i P_i,
\end{equation}
where in general the orthogonal projections $P_i$ may be of higher rank than one.
By plugging it into (\ref{App:5}) one arrives at
\begin{equation}
    \sum_{i=0}^{d-1}\mathrm{Re}\left(\omega^{(i+j)k}\langle b_j|P_i|b_j\rangle\right)\geq 1-2d\sqrt{2\varepsilon},
\end{equation}
which after summing over $k=0,\ldots,d-1$ gives 
\begin{equation}\label{App:6}
    \langle b_j|P_{-j}|b_j\rangle\geq 1-2d\sqrt{2\varepsilon}.
\end{equation}
This inequality means that each vector $\ket{b_j}$ is $\varepsilon$-close to 
the subspace corresponding to the $\omega^{-j}$ outcome of $B_2$. Let us then 
introduce vectors $\ket{v_j}=P_{-j}\ket{b_j}$ and their normalized versions
$\ket{\overline{v}_j}=\ket{v_j}/\|\ket{v_j}\|$. From the fact that the projections $P_j$ are pairwise orthogonal we infer that $\ket{\overline{v}_j}$ are orthogonal as well. Moreover, Eq. (\ref{App:6}) implies that $\|\ket{v}_j\|\geq 1-2d\sqrt{2\varepsilon}$, and so the norm of each $\ket{v_j}$ is close to unity.

Now, due to the fact that by the very construction each $\ket{\overline{v}_j}$
belongs to the subspace corresponding to $P_{-j}$, the observable $B_2$ can be written as
\begin{equation}
    B_2=\sum_{i=0}^{d-1}\omega^{-i}\proj{\overline{v}_i}\oplus B_2',
\end{equation}
where $B_2'$ is an operator whose support is orthogonal to 
all $\ket{v_i}$. Then, due to the fact that $\ket{v_i}$ are orthogonal
one can find a unitary operation such that $U_B\ket{\overline{v}_i}=\ket{d-i}$. This operation brings $B_2$ to
\begin{equation}\label{App:65}
    U_B B_2 U_B^{\dagger}=\sum_{i=0}^{d-1}\omega^i\proj{i}\oplus B_2''.
\end{equation}
Denoting then $\ket{b_i'}=U_B\ket{b_i}$, where $\ket{b_i}$ are vectors appearing in the decomposition of $\ket{\psi}$ in Eq. (\ref{stateRob}), we deduce from
Eq. (\ref{App:6}) that
\begin{equation}\label{App:7}
    \langle b_i'|i\rangle\geq 1-2d\sqrt{2\varepsilon}.
\end{equation}

Given conditions (\ref{App:7}) we can now move on to proving 
our main statements (\ref{Rob1}) and (\ref{Rob2}). To this end, 
we notice that for the unitary operation $U_B$ introduced above, one has
\begin{eqnarray}\label{App:75}
\left\|(\mathbbm{1}\otimes U_BB_i^k)\ket{\psi}-[\mathbbm{1}\otimes (B'_i)^{k}]\ket{\phi^+_d}\right\|&=&
\left\|(A_i^{k}\otimes U_B B_i^{k})\ket{\psi}-[A_i^{k}\otimes (B'_i)^{k}]\ket{\phi^+_d}\right\|\nonumber\\
&=&\left\|(A_i^{k}\otimes \widetilde{B}_i^{k})\ket{\widetilde{\psi}}-\ket{\widetilde{\psi}}+\ket{\widetilde{\psi}}-\ket{\phi^+_d}\right\|,
\end{eqnarray}
where the first equation is a consequence of the fact that 
$A_i^k$ are unitary and that the vector norm is unitarily invariant, whereas the second equation has been obtained by adding and subtracting $\ket{\widetilde{\psi}}$ with the latter defined as $\ket{\widetilde{\psi}}=(\mathbbm{1}\otimes U_B)\ket{\psi}$ and 
$\widetilde{B}_i=U_B B_iU_B^{\dagger}$; in particular $\widetilde{B}_2$ is given in Eq. (\ref{App:65}). 

Using then the triangle inequality we arrive at
\begin{eqnarray}\label{App:8}
\left\|(A_i^{k}\otimes \widetilde{B}_i^{k})\ket{\widetilde{\psi}}-\ket{\widetilde{\psi}}+\ket{\widetilde{\psi}}-\ket{\phi^+_d}\right\|\leq\left\|(A_i^{k}\otimes \widetilde{B}_i^{k})\ket{\widetilde{\psi}}-\ket{\widetilde{\psi}}\right\|+\left\|\ket{\widetilde{\psi}}-\ket{\phi^+_d}\right\|
\end{eqnarray}
Let us now consider the first term of the right-hand side of the above inequality
and write is as
\begin{eqnarray}\label{App:85}
    \left\|(A_i^{k}\otimes \widetilde{B}_i^{k})\ket{\widetilde{\psi}}-\ket{\widetilde{\psi}}\right\|&=&\left\|(A_i^{k}\otimes B_i^{k})\ket{\psi}-\ket{\psi}\right\|\nonumber\\
    &=&\left\{2\left[1-\mathrm{Re}\left(\langle \psi|A_i^{k}\otimes B_i^{k}|\psi\rangle\right)\right]\right\}^{1/2}\nonumber\\
    &\leq &\sqrt{2\varepsilon},
\end{eqnarray}
where to get the inequality we have used conditions (\ref{App:1}). 

Let us then consider the second term appearing on the right-hand side of 
Eq. (\ref{App:8}) and rewrite it as
\begin{eqnarray}\label{App:10}
    \left\|\ket{\widetilde{\psi}}-\ket{\phi^+_d}\right\|&=&\left\{2\left[1-\mathrm{Re}(\langle\widetilde{\psi}|\phi_{d}^+\rangle)\right]\right\}^{1/2}\nonumber\\
    &=&\left\{2\left[1-\frac{1}{\sqrt{d}}\sum_{i}\alpha_i\,\mathrm{Re}
    (\langle\overline{b}_i|i\rangle)\right]\right\}^{1/2}\nonumber\\
    &\leq &\left\{2\left[1-\frac{1}{\sqrt{d}}(1-2d\sqrt{2\varepsilon})\sum_{i}\alpha_i\right]\right\}^{1/2}\nonumber\\
    &\leq & \left\{2\left[1-(1-2d\sqrt{2\varepsilon})\sqrt{1-
    \varepsilon}\right]\right\}^{1/2}\nonumber\\
    &\leq & 2\sqrt{d}\sqrt[4]{2\varepsilon},
\end{eqnarray}
where to obtain the first inequality we have used Eqs. (\ref{App:7})
and (\ref{App:2}), respectively, whereas the third inequality stems from 
two simple inequalities: $1-\sqrt{1-\varepsilon}\leq \varepsilon$ and $\varepsilon\leq \sqrt{2\varepsilon}$. By applying Eqs. (\ref{App:8}), (\ref{App:85}) and (\ref{App:10}) into (\ref{App:75}) we obtain (\ref{Rob1}).


Let us eventually prove (\ref{Rob2}). For this purpose, we first observe that
\begin{equation}\label{App:9}
  \left\|\widetilde{B}_i^{k}-(B'_i)^{k}\right\|_2^2=d\left\|\left[\widetilde{B}_i^{k}-(B'_i)^{k}\right]\ket{\phi_d^{+}}\right\|^2.
\end{equation}
Then, the triangle inequality gives us
\begin{eqnarray}
    \left\|\widetilde{B}_i^{k}\ket{\widetilde{\psi}}-(B'_i)^{k}\ket{\phi_d^+}\right\|&=&
    \left\|\widetilde{B}_i^{k}\ket{\widetilde{\psi}}-(B'_i)^{k}\ket{\phi_d^+}+(B'_i)^k\ket{\phi_d^+}-\widetilde{B}_i^k\ket{\phi_d^+}\right\|\nonumber\\
    &\geq&\left\|\left[\widetilde{B}_i^{k}-(B'_i)^k\right]\ket{\phi_d^+}\right\|-\left\|\widetilde{B}_i^{k}(\ket{\widetilde{\psi}}-\ket{\phi_d^+})\right\|,
\end{eqnarray}
which by taking into account the identity (\ref{App:9}) and the unitary invariance of the vector norm, leads us to
\begin{equation}
    \left\|\widetilde{B}_i^{k}-(B'_i)^{k}\right\|_2^2\leq  \sqrt{d}\left(\left\|\widetilde{B}_i^{k}\ket{\widetilde{\psi}}-(B'_i)^{k}\ket{\phi_d^+}\right\|+\left\|\ket{\widetilde{\psi}}-\ket{\phi_d^+}\right\|\right).
\end{equation}
Using inequalities (\ref{Rob1}) and (\ref{App:10}) we obtain 
(\ref{Rob2}).
\end{proof}




\end{document}