Stationarity, discrete-translation operator, and the power spectral density matrix
Given the definition of the correlation matrix $\mathbf{R}_{\mathbf{x}}$ here, I am assuming that $\mathsf{E}[\mathbf{x}] = \mathbf{0}$. I do this because the correlation matrix is usually defined as $\mathsf{E}[(\mathbf{x} - \mathsf{E}[\mathbf{x}])(\mathbf{x} - \mathsf{E}[\mathbf{x}])^{\dagger}]$, where $\dagger$ indicates complex conjugate tranpose.
 
Note that since $\mathbf{R}_{\mathbf{x}}$ is a correlation matrix, it is Hermitian positive semi-definite, which means that it is unitarily diagonalizable and has all non-negative eigenvalues:
\begin{equation}
\mathbf{R}_{\mathbf{x}} = \mathbf{U}\mathbf{\Lambda}\mathbf{U}^{\dagger},
\end{equation}
where $\mathbf{U}$ is a unitary matrix and $\mathbf{\Lambda}$ is a diagonal matrix with non-negative numbers on its diagonal.

Let's consider the translation operator $\mathbf{T}$:
\begin{equation}
\mathbf{T} = \mathbf{\Phi}\mathbf{P}\mathbf{\Phi}^{\dagger},
\end{equation}
where $\mathbf{\Phi}^{\dagger}$ is the DFT matrix and $\mathbf{P}$ is diagonal. As I mentioned in a comment, $\mathbf{T}$ should be unitary, and that requires that the diagonal matrix $\mathbf{P}$ have complex exponentials on its diagonal. In fact, to make $\mathbf{\Phi}\mathbf{P}\mathbf{\Phi}^{\dagger}$ a translation operator, $\mathbf{P}$ should be
\begin{equation}
\mathbf{P} = \textrm{diag}(1,e^{2\pi i/M},\ldots,e^{2\pi i(M-1)/M}).
\end{equation}


I do not think this will have an impact on the proof, but it is important for applications.

Now let's assume that $\mathbf{R}_{\mathbf{T}\mathbf{x}} = \mathbf{R}_{\mathbf{x}}$:
\begin{eqnarray}
\mathsf{E}[(\mathbf{T}\mathbf{x})(\mathbf{T}\mathbf{x})^{\dagger}] &amp;=&amp; \mathbf{R}_{\mathbf{x}}\\
\mathbf{T}\mathsf{E}[\mathbf{x}\mathbf{x}^{\dagger}]\mathbf{T}^{\dagger} &amp;=&amp; \mathbf{R}_{\mathbf{x}}\\
\mathbf{T}\mathbf{R}_{\mathbf{x}}\mathbf{T}^{\dagger} &amp;=&amp; \mathbf{R}_{\mathbf{x}}\\
\mathbf{T}\mathbf{R}_{\mathbf{x}} &amp;=&amp; \mathbf{R}_{\mathbf{x}}\mathbf{T}
\end{eqnarray}
Since


$\mathbf{R}_{\mathbf{x}}$ and $\mathbf{T}$ are both diagonalizable
and they commute with each other,


they are simultaneously diagonalizable. This means that there is a single matrix that diagonalizes both.

We have been told from the beginning that $\mathbf{\Phi}$ diagonalizes $\mathbf{T}$, so now we know that $\mathbf{\Phi}$ diagonalizes $\mathbf{R}_{\mathbf{x}}$, too. This means that the unitary matrix $\mathbf{U}$ that diagonalizes $\mathbf{R}_{\mathbf{x}}$ must be $\mathbf{\Phi}$:
\begin{equation}
\mathbf{R}_{\mathbf{x}} = \mathbf{\Phi}\mathbf{\Lambda}\mathbf{\Phi}^{\dagger}.
\end{equation}
We have already established that $\mathbf{\Lambda}$ is diagonal with non-negative real numbers on its diagonal. The $\mathbf{S}_{\mathbf{x}}$ that we have sought, is $\mathbf{\Lambda}$.