Before delving into the theory of the covariance matrix criterion, let us first discus the definition and properties of the covariant matrices.

### Definition and properties

**Def. 5** *Given a quantum state $\rho$ and a complete set of orthonormal observables $\{M_k\}$ on a $d$-dimensional Hilbert space, the $d^2 \times d^2$ covariant matrix $\gamma = \gamma(\rho,\{M_k\})$ and symmetrized covariant matrix $\gamma^S = \gamma^S(\rho,\{M_k\})$ are defined as $$\begin{aligned}
\gamma_{ij} &=& \langle M_i M_j \rangle - \langle M_i \rangle \langle M_j \rangle, \\
\gamma^S_{ij} &=& \frac{\langle M_i M_j \rangle + \langle M_j M_i \rangle}{2} - \langle M_i \rangle \langle M_j \rangle,\end{aligned}$$ where $\langle M \rangle = \mbox{tr}(\rho M)$ defines the expectation of the operator $M$ with respect to the state $\rho$.*

The complete set of orthonormal observables $\{M_k\}$ has to satisfy the Hilbert-Schmidt orthonormality condition $\mbox{tr}(M_i M_j) = \delta_{ij},$ $\forall i,j = 1,2,...,d$. One example for such a set of observables for the case of single qubit in terms of the Pauli matrices, can be given by $$M_1 = \frac{I}{\sqrt{2}}, \ M_2 = \frac{\sigma_x}{\sqrt{2}}, \ M_3 = \frac{\sigma_y}{\sqrt{2}},
\ M_4 = \frac{\sigma_z}{\sqrt{2}}.$$ In general, for the $d$-dimensional case, one can consider following matrices to form the complete set of orthonormal observables: $$\begin{aligned}
X_i 