Local symplectic transformations on Gaussian states According to Eqn. $(18)$ of this paper, given a two mode Gaussian state with the $4 \times 4$ covariance matrix $\sigma$, it is possible to find a symplectic matrix $S = S_1 \oplus S_2$, where $S_{1, 2}$ are local $2 \times 2$ symplectic matrices, such that:
$$\sigma' = S \sigma S^T = \begin{pmatrix}a & 0 & c_1 & 0
\\
 0 & a & 0 & c_2
\\ c_1 & 0 & b & 0 \\
 0 & c_2 & 0 & b
 \end{pmatrix} $$
Where $a, b, c_{1, 2}$ are real. How to show this? I can see how the diagonal blocks can become diagonal, from the fact that the covariance matrices can be diagonalized using symplectic matrices. But how come the off diagonal terms are diagonalized? Is there any reference which derives this result?
If $$S = \begin{bmatrix}S_1 & 0 \\ 0 & S_2 \end{bmatrix}$$ and $$\sigma = \begin{bmatrix}\sigma_A & \sigma_{AB} \\ \sigma_{AB}^T & \sigma_B \end{bmatrix},$$ then after the transformation
$$\sigma' = 
\begin{bmatrix}
S_1 \sigma_A S_1^T & S_1 \sigma_{AB} S_2^T \\
S_2 \sigma_{AB}^T S_1^T & S_2 \sigma_B S_2^T
\end{bmatrix}
$$
According to Williamson's theorem, we can choose $S_{1, 2}$ to diagonalize $\sigma_{A, B}$, so the diagonal blocks can be made diagonal. What about the off diagonal terms?
Moreover, can this be extended to larger systems, say $4$ mode systems with $8 \times 8$ covariance matrices?
 A: The trick here is that two-dimensional orthogonal matrices (rotations) are also symplectic matrices. They correspond to phase shifts of phase space, as explained more carefully here [1].
\begin{equation} 
R(\phi) = \begin{pmatrix} \cos(\phi) & -\sin(\phi) \\
\sin(\phi) & \cos(\phi) \end{pmatrix}
\end{equation}
So, after applying $S_1 \oplus S_2$ to bring $\sigma$ to the form
\begin{equation}
\sigma = \begin{pmatrix} a 1_{2x2} & C \\ C^T & b 1_{2x2} \end{pmatrix}
\end{equation}
you can then decompose $C$ using the singular value decomposition as $C = R(\phi_1) D R( - \phi_2)$, where $D = \begin{pmatrix} c_1 & 0 \\ 0 & c_2\end{pmatrix}$ and then apply the symplectic transformation $R(\phi_1) \oplus R(\phi_2)$ to bring the covariance matrix to your desired form. I do not think that this works exactly like this for larger systems since the orthogonal matrices are contained in the symplectic matrices only for the special case of two dimensions.
[1] https://arxiv.org/pdf/1401.4679.pdf
A: This is another answer, less mathematical than the one by @saidthemonkey, but basically equivalent. You can see it as a more physical point of view of the process to make the covariance matrix in the $\sigma'$ form. This process can be decomposed in the following steps

*

*Diagonalize the diagonal blocks using orthogonal matrices = rotation matrices = dephasing;

*Equalize the two values of each of the block through a diagonal symplectic matrix $\begin{bmatrix} s&0\\0& \frac1s\end{bmatrix}$ = a squeezing tranformation;

*Since both diagonal blocks are now $\propto I$, they are invariant under rotations. So you can use further rotations to diagonalize the off-diagonal blocks

