\documentclass[11pt]{article}
\usepackage{latexsym}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{amsthm}
\usepackage{hyperref}
\usepackage{braket}
\usepackage{xcolor}
%\usepackage{epsfig}
%\usepackage{psfig}

\newcommand{\handout}[5]{
  \noindent
  \begin{center}
  \framebox{
    \vbox{
      \hbox to 5.78in { {\bf PHYS 7895: Gaussian Quantum Information } \hfill #2 }
      \vspace{4mm}
      \hbox to 5.78in { {\Large \hfill #5  \hfill} }
      \vspace{2mm}
      \hbox to 5.78in { {\em #3 \hfill #4} }
    }
  }
  \end{center}
  \vspace*{4mm}
}

\newcommand{\lecture}[4]{\handout{#1}{#2}{#3}{Scribe: #4}{Lecture #1}}

%\newtheorem{theorem}{Theorem}
%\newtheorem{corollary}[theorem]{Corollary}
%\newtheorem{lemma}[theorem]{Lemma}
%\newtheorem{observation}[theorem]{Observation}
%\newtheorem{proposition}[theorem]{Proposition}
%\newtheorem{definition}[theorem]{Definition}
%\newtheorem{claim}[theorem]{Claim}
%\newtheorem{fact}[theorem]{Fact}
%\newtheorem{assumption}[theorem]{Assumption}
%
%\usepackage{amsthm}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{prop}[theorem]{Proposition}
\newtheorem{cor}[theorem]{Corollary}
\theoremstyle{definition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{remark}[theorem]{Remark}
\newtheorem{example}[theorem]{Example}

% 1-inch margins, from fullpage.sty by H.Partl, Version 2, Dec. 15, 1988.
\topmargin 0pt
\advance \topmargin by -\headheight
\advance \topmargin by -\headsep
\textheight 8.9in
\oddsidemargin 0pt
\evensidemargin \oddsidemargin
\marginparwidth 0.5in
\textwidth 6.5in

\parindent 0in
\parskip 1.5ex


\newcommand{\brac}[1]{\lbrace #1 \rbrace}
\newcommand{\mi}{{\min}}
\newcommand{\ma}{{\max}}
\def\U{\mathrm{U}}
\def\A{\mathcal{A}}
\def\B{\mathcal{B}}
\def\C{\mathcal{C}}
\def\D{\mathcal{D}}
\def\E{\mathcal{E}}
\def\F{\mathcal{F}}
\def\G{\mathcal{G}}
\def\H{\mathcal{H}}
\def\I{\mathcal{I}}
\def\L{\mathcal{L}}
\def\M{\mathcal{M}}
\def\N{\mathcal{N}}
\def\O{\mathcal{O}}
\def\P{\mathcal{P}}
\def\Q{\mathcal{Q}}
\def\T{\mathcal{T}}
\def\R{\mathbb{R}}
\def\J{\mathcal{J}}
\def\S{\mathcal{S}}
\def\K{\mathcal{K}}
\def\X{\mathcal{X}}
\newcommand{\Tr}{\operatorname{Tr}}
\newcommand{\Det}{{\operatorname{Det}}}
\newcommand{\unit}{1\!\!1}
\newcommand\rit[1]{\textcolor{red}{\textit{#1}}}
\newcommand\bit[1]{\textcolor{blue}{\textbf{#1}}}
\newcommand\pit[1]{\textcolor{pi}{\textbf{#1}}}
\makeatletter
\let\@fnsymbol\@arabic
\makeatother

%\renewcommand{\baselinestretch}{1.25}

\begin{document}

\lecture{9}{}{Lecturer: Mark M.\ Wilde}{Vishal Katariya}

This document is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 Unported License.

%\tableofcontents

\section{Overview}

In the previous lecture, we defined and studied faithful Gaussian states as thermal states of quadratic Hamiltonians.

In this lecture, we continue the analysis of faithful states as thermal states of quadratic Hamiltonians, and via the Williamson Theorem, we show and prove the form of a general Gaussian state.
 
%We point readers to \cite{AS17} for background on topics covered in this lecture. 

\section{Recap}

In the previous lecture, we defined a faithful Gaussian state to be a thermal state
\begin{equation}
	\frac{e^{-  \hat{H}}}{\Tr[e^{-  \hat{H}}]} 
\end{equation}
of a quadratic Hamiltonian
\begin{equation}
	\hat{H} = \frac{1}{2} (\hat{r} - \overline{r})^T H (\hat{r} - \overline{r})
\end{equation}
where $\overline{r} \in \R^{2n}$ and $H$ is a $2n \times 2n$ positive definite real matrix.

We showed how to build up a faithful Gaussian state with Hamiltonian matrix 
\begin{align}
H' &  = S^T H S, \label{eq:H-prime}\\
H & = \bigoplus_{j=1}^n \lambda_j \left( \begin{matrix} 1&0\\ 0&1 \end{matrix} \right)
\end{align}
 with $\lambda_j > 0 \ \forall j \in \{1, \ldots, n\}$, $S$ symplectic, and Hamiltonian operator
\begin{align}
	\hat{H}' &= \frac{1}{2} (\hat{r} - \overline{r})^T S^T H S (\hat{r} - \overline{r}) \\
			 &= \frac{1}{2} (\hat{r} - \overline{r})^T H' (\hat{r} - \overline{r}).
\end{align}
This yields a quantum Gaussian state with the following density operator:
\begin{equation}
	\rho_G = \frac{e^{-\frac{1}{2}(\hat{r} - \overline{r})^T H' (\hat{r} - \overline{r})}}{\sqrt{\Det\left( \frac{\sigma' + i \Omega}{2} \right)}} \label{eq:general-rho}
\end{equation}
where the mean vector of $\rho_G$ is $\overline{r}$ and the covariance matrix is 
\begin{equation}
	\sigma' = S^{-1} \sigma S^{-T} \quad \text{ with } \quad \sigma = \bigoplus_{j=1}^n \coth\left(\frac{\lambda_j}{2}\right) \left( \begin{matrix} 1&0\\ 0&1 \end{matrix} \right). \label{eq:cov-mat-from-bef}
\end{equation}

Also, we noted that $\rho_G$ can be written as
\begin{equation}
	\rho_G = \frac{\hat{D}_{- \overline{r} } \hat{S} e^{-\frac{1}{2}\hat{r} ^T H' \hat{r}} \hat{S}^{\dag} \hat{D}_{ \overline{r} } }{\sqrt{\Det\left( \frac{\sigma' + i \Omega}{2} \right)}} \label{eq:rho-symp}
\end{equation}
where $S = e^{\frac{i}{2} \hat{r}^T (\Omega^T \ln S )\hat{r} }  $ and $\hat{D}_{ \overline{r} } = \exp\left( i \overline{r}^T \Omega \hat{r} \right)$.\footnote{Note that for some $S$, we may need two quadratic evolutions, not one.}

Since $H$ is diagonal, 
\begin{equation}
	\frac{  e^{-\frac{1}{2}\hat{r} ^T H \hat{r}} }{ \sqrt{\Det\left( \frac{\sigma' + i \Omega}{2} \right)} } = \bigotimes_{j=1}^n \frac{ e^{-\lambda_j(\hat{n} + \frac{1}{2})} }{z(\lambda_j)} 
\end{equation}
with $z(\lambda_j) = \left[e^{\lambda_j/ 2} - e^{-\lambda_j/2}\right]^{-1}$.

Note that $\frac{ e^{-\lambda(\hat{n} + \frac{1}{2})} }{z(\lambda)}$ is typically called the bosonic thermal state, and has mean photon number $\langle\hat{n}\rangle = \frac{1}{2} \langle \hat{x}^2 + \hat{p}^2 - 1 \rangle = \coth\left(\frac{\lambda}{2}\right) - \frac{1}{2}$.

\section{Form of a general Gaussian state}

Now we will prove, perhaps surprisingly, that the state given in \eqref{eq:general-rho} is the most general form that a faithful Gaussian state can take.

\begin{theorem}
	A Gaussian state given in the form
	\begin{equation}
	\rho_G = \frac{e^{-\frac{1}{2}(\hat{r} - \overline{r})^T H' (\hat{r} - \overline{r})}}{\sqrt{\Det\left( \frac{\sigma' + i \Omega}{2} \right)}},
	\end{equation}
	with $H'$ as in \eqref{eq:H-prime},
	 is the most general form for any faithful Gaussian state.
\end{theorem}

%\begin{lemma}
%It suffices to focus on Hamiltonians of the form 
%\begin{equation}
%	H' = \frac{1}{2} (\hat{r} - \overline{r})^TH(\hat{r} - \overline{r})
%\end{equation}
%when defining faithful Gaussian states.
%\end{lemma}

\begin{proof}
Suppose that the faithful Gaussian state is given by
\begin{equation}
	\rho_G = \frac{ e^{ - \frac{1}{2} (\hat{r} - \overline{r})^T H (\hat{r} - \overline{r})} }{\Tr\left[ e^{ - \frac{1}{2} (\hat{r} - \overline{r})^T H (\hat{r} - \overline{r})} \right]}
\end{equation}
where $\overline{r} \in \R^{2n}$ and $H$ is a $2n \times 2n$ real positive definite matrix.

Then, by the Williamson Theorem given below, there exists symplectic $S$ such that
\begin{equation}
	H = S^{T} H_{\text{diag}} S
\end{equation}
where $H_{\text{diag}} = \bigoplus_{j=1}^n \lambda_j \left( \begin{matrix} 1&0 \\ 0&1 \end{matrix} \right)$ for $\lambda_j > 0 \quad \forall j \in \{1,\ldots,n\}$
so that we can write $\rho_G$ as
\begin{equation}
\frac{ e^{ - \frac{1}{2} (\hat{r} - \overline{r})^T S^{T} H_{\text{diag}} S (\hat{r} - \overline{r})} }{\Tr\left[ e^{ - \frac{1}{2} (\hat{r} - \overline{r})^T S^{T} H_{\text{diag}} S (\hat{r} - \overline{r})} \right]}
\end{equation}
Then by manipulations we have already conducted, there exists a unitary $\hat{S}$ (generated by a quadratic Hamiltonian) and a displacement operator $\hat{D}_{ \overline{r}} = \exp(i \overline{r}^T \Omega \hat{r})$ such that $\rho_G$ is given as
\begin{equation}
\propto \hat{D}_{- \overline{r} } \hat{S} e^{-\frac{1}{2}\hat{r} ^T H' \hat{r}} \hat{S}^{\dag} \hat{D}_{ \overline{r} } 
\end{equation}
and the normalization is given by
\begin{equation}
\sqrt{\Det\left( \frac{\sigma' + i \Omega}{2} \right)},
\end{equation}
with $\sigma'$ defined in \eqref{eq:cov-mat-from-bef}.
%
%
%Therefore it suffices to focus on Hamiltonians of the form 
%\begin{equation}
%	H' = \frac{1}{2} (\hat{r} - \overline{r})^TH(\hat{r} - \overline{r})
%\end{equation}
%when defining faithful Gaussian states.
\end{proof}

%Consider that 
%\begin{align}
%	\hat{H}' &= \frac{1}{2} \sum_{jk} ( \hat{r}_j - \overline{r}) H_{jk} ( \hat{r}_j - \overline{r}_k) \\
%			 &= \frac{1}{2} \sum_{jk} H_{jk} ( \hat{r}_j - \overline{r}_k)( \hat{r}_k - \overline{r}_k ). \label{eq:ham-comp}
%\end{align}
%
%Last time we showed that 
%\begin{equation}
%	\hat{D}_{ - \overline{r} } \hat{r} \hat{D}_{\overline{r}} = \hat{r} - \overline{r}
%\end{equation}
%where $\hat{D}_{ \overline{r}} = \exp(i \overline{r}^T \Omega \hat{r})$.
%
%So this means that upon substituting the above into \eqref{eq:ham-comp}, we get
%\begin{align}
%	\hat{H}' &= \frac{1}{2} \sum_{jk} H_{jk} \left( \hat{D}_{ \overline{-r}} \hat{r}_j \hat{D}_{ \overline{r}} \right) \left( \hat{D}_{- \overline{r}} \hat{r}_k \hat{D}_{ \overline{r}} \right) \\
%			 &= \frac{1}{2} \sum_{jk} H_{jk} \hat{D}_{ - \overline{r}} \hat{r}_j \hat{r}_k \hat{D}_{ \overline{r}} \\
%			 &= \hat{D}_{ - \overline{r}} \left[ \frac{1}{2} \hat{r}^T H \hat{r} \right] \hat{D}_{ \overline{r}}.
%\end{align}
%This means that we can realize the shift in quadrature operators of Hamiltonian $\hat{H}$ via a unitary displacement operator.

\subsection{Williamson Theorem}

\begin{theorem}[Williamson]
	Given a $2n \times 2n$ positive definite real matrix $M$, there exists a symplectic transformation $S$ such that
	\begin{equation}
		S M S^T = D ,
	\end{equation}
	with
	\begin{equation}
	D = \bigoplus_{j=1}^n d_j \left( \begin{matrix} 1&0\\ 0&1 \end{matrix} \right)
	\end{equation}
	 and $  d_j > 0 \ \forall j \in \{1,\ldots,n\}$. The set $\{d_j\}_{j=1}^n$ is the set of symplectic eigenvalues of $M$. 
\end{theorem}

Before proving the Williamson Theorem, we recall a standard lemma about the decomposition of real antisymmetric matrices.

\begin{lemma} \label{lem:antisymm}
	Let $A$ be a real, full-rank, antisymmetric $2n \times 2n$ matrix (i.e., $A=-A^T$). Then there exists a real orthogonal $2n \times 2n$ matrix $O$ such that
	\begin{equation}
		O A O^T = \bigoplus_{j=1}^n c_j \Omega_1,
	\end{equation}
	where $\Omega_1 = \left( \begin{matrix} 0&1 \\ -1&0 \end{matrix} \right)$ and $c_j > 0$.
\end{lemma}

\begin{proof}
	Since $A$ is antisymmetric and full rank, it follows that the matrix $A^2$ is symmetric and negative definite because
	\begin{equation}
	A^2 = A A = - A^T A < 0
	\end{equation}
	since $A^T A$ is positive definite for any full rank $A$.
	Thus there exists an orthogonal transformation $O'$ such that $O' A^2 O'^T = B$ with $B$ diagonal and having strictly negative entries.

	Let $\ket{\psi}$ be some eigenvector of $A^2$ with eigenvalue $b_1 < 0$. 

	Then
	\begin{equation}
	\left\|A \ket{\psi}\right\|^2 = \bra{\psi} A^T A \ket{\psi} = - \bra{\psi} A^2 \ket{\psi} = -b_1 = |b_1|.
	\end{equation}

	So $\ket{\psi'} = \frac{A \ket{\psi}}{\sqrt{|b_1|}}$ is normalized and  orthogonal to $\ket{\psi}$ because
	\begin{align}
		\braket{\psi | \psi'} &= \frac{\braket{\psi |A \psi}}{\sqrt{|b_1|}} =
		\frac{\left( \braket{\psi | A \psi}\right)^T}{\sqrt{|b_1|}} \\
							 &= \frac{\braket{\psi | A^T \psi}}{\sqrt{|b_1|}}
							 = - \frac{ \braket{\psi | A \psi} }{\sqrt{|b_1|}} = 0.
	\end{align}
	In the above, we showed that $\braket{\psi | \psi'} = - \braket{\psi | \psi'}$ and thus  $\braket{\psi | \psi'} = 0$.

	Suppose now that $\ket{\phi}$ is in the subspace orthogonal to $\text{span} \{ \ket{\psi}, \ket{\psi'} \}$. This implies that 
	\begin{align}
		\braket{\phi | A\psi} &= \braket{\phi | \psi'} \sqrt{|b_1|} = 0, \\
		\braket{\phi | A\psi'} &= \frac{\braket{\phi | A^2 \psi}}{\sqrt{|b_1|}} = \frac{\braket{\phi | \psi}b_1}{\sqrt{|b_1|}} \\
							  &= - \braket{\phi | \psi}\sqrt{|b_1|} \\
							  &= 0.
	\end{align}
	Furthermore, due to the antisymmetry of $A$, $\braket{\psi | A \psi} = 0 = \braket{\psi' | A \psi'}$.

	Also,
	\begin{align}
		\braket{\psi | A \psi'} &= \frac{\braket{\psi | A^2 \psi}}{\sqrt{|b_1|}} = \frac{b_1 \braket{\psi | \psi}}{\sqrt{|b_1|}} 
							   = -\sqrt{|b_1|}, \\
		\braket{\psi' | A \psi} &= \sqrt{|b_1|},
	\end{align}
		where the second statement above is due to the antisymmetry of $A$.
		
		Now define the orthogonal matrix $O_1$ as
		\begin{equation}
		O_1 = \begin{bmatrix}
		    \ket{\psi'} & \ket{\psi} & \ket{v_1} & \ldots & \ket{v_{2n-2}}
		\end{bmatrix},
		\end{equation}
		where $\ket{v_1}, \ldots, \ket{v_{2n-2}}$ is a set of orthonormal vectors orthogonal to $\ket{\psi'}$ and $ \ket{\psi}$.

	Putting everything above together, we conclude that
	\begin{equation}
		O_1^T A O_1 = \begin{pmatrix} \sqrt{|b_1|} \Omega_1&0 \\ 0&A' \end{pmatrix} 
	\end{equation}
	and we see that this step gives the first step of the decomposition, after setting $c_1 = \sqrt{|b_1|}$. The matrix $A'$ is antisymmetric, and so this procedure can be repeated exhaustively to complete the decomposition. 

	%\textit{Note}: In the above, we took $b_1 > 0$. However, that can be done with no loss in generality because $\begin{pmatrix} 0 & -b_1 \\ b_1&0 \end{pmatrix}$  is related to $\begin{pmatrix} 0 & b_1 \\ -b_1&0 \end{pmatrix}$ by an orthogonal similarity transformation described by the matrix $\begin{pmatrix} 0 & 1 \\ 1&0 \end{pmatrix}$.
\end{proof}

\begin{proof} \textbf{Proof of Williamson Theorem}

	Consider the matrix $M^{-\frac{1}{2}} \Omega M^{-\frac{1}{2}}$. It is real and $2n \times 2n$. 
	
	Due to the symmetry of $M$ and antisymmetry of $\Omega$, it follows that $M^{-\frac{1}{2}} \Omega M^{-\frac{1}{2}}$ is antisymmetric.

	Additionally, since both $M$ and $\Omega$ are full rank, $M^{-\frac{1}{2}} \Omega M^{-\frac{1}{2}}$ is also full rank.

	Thus by invoking Lemma~\ref{lem:antisymm} above, there exists a real orthogonal transform $O$ such that
	\begin{equation}
		O M^{-\frac{1}{2}} \Omega M^{-\frac{1}{2}} O^T = \bigoplus_{j=1}^{n} d_j^{-1} \Omega_1 \quad \text{with} \quad d_j > 0.
	\end{equation}

   Define $\underline{D} = \begin{bmatrix} d_1 &  &  \\  & \ddots & \\ & & d_n \end{bmatrix}$, an $n\times n$ matrix, and  set $D = \underline{D} \otimes I_2.$

   Then we have $\bigoplus_{j=1}^n \frac{1}{d_j} \Omega_1 = \underline{D}^{-1} \otimes \Omega_1$. Using this expression, we find that
   \begin{align}
	   D^{\frac{1}{2}} \left[\bigoplus_{j=1}^n \frac{1}{d_j} \Omega_1 \right] D^{\frac{1}{2}} 
	   & = \left( \underline{D}^{\frac{1}{2}} \otimes I_2\right) \left( \underline{D}^{-1} \otimes \Omega_1 \right) \left( \underline{D}^{\frac{1}{2}} \otimes I_2\right) \\
	   & = I_n \otimes \Omega_1 = \Omega.
   \end{align}
   \begin{equation}
   		\implies D^{\frac{1}{2}} O M^{- \frac{1}{2}} \Omega M^{ - \frac{1}{2}} O^T D^{ \frac{1}{2}} = \Omega.	
		\label{eq:williamson-symp-constr}
   \end{equation}
     Now set
     \begin{equation}
     S \equiv D^{ \frac{1}{2}} O M^{- \frac{1}{2}},
     \end{equation}
     and we conclude from \eqref{eq:williamson-symp-constr} that $S \Omega S^T  = \Omega$, so that $S$ is symplectic.
  
   Also,
   \begin{align}
	   SMS^T &= \left( D^{ \frac{1}{2}} O M^{- \frac{1}{2}} \right) M \left( D^{ \frac{1}{2}} O M^{ - \frac{1}{2}} \right)^T \\
	   &= D^{ \frac{1}{2}} O M^{- \frac{1}{2}} M M^{- \frac{1}{2}} O^T D^{ \frac{1}{2}} \\
	   &= D^{ \frac{1}{2}} O O^T D^{ \frac{1}{2}} \\
	   &= D^{ \frac{1}{2}} D^{ \frac{1}{2}} \\
	   &= D.
   \end{align}
   The statement of the Williamson Theorem is that for a positive definite and real $2n \times 2n$ matrix $M$, there exists a $2n \times 2n$ symplectic matrix $S$ such that $SMS^T = D$. We have proved by construction of $S$ that such a transform exists. Thus, we have completed the proof of the Williamson Theorem.   
\end{proof}

We can now apply the Williamson Theorem to the Hamiltonian matrix for a faithful Gaussian state.

We recall the form of a Gaussian state
\begin{align}
	\frac{e^{-\hat{H}}}{\Tr\big[e^{-\hat{H}}\big]} &= \frac{e^{- \frac{1}{2} ( \hat{r} - \overline{r})^T H ( \hat{r} - \overline{r})}}{\Tr[e^{-\hat{H}}]} \\
												   &= \frac{\hat{D}_{ - \overline{r}} e^{- \frac{1}{2} \hat{r}^T H \hat{r}} \hat{D}_{ \overline{r}} }{\Tr\big[ e^{- \frac{1}{2} \hat{r}^T H \hat{r}} \big]}.
\end{align}

   We use the symplectic diagonalization of $H$ as $H = S_H^T \left( \Lambda \otimes I_2 \right) S_H$ where $\Lambda = \begin{bmatrix} \lambda_1 &  &  \\  & \ddots & \\ & & \lambda_j \end{bmatrix}$, and $S_H$ is the transposed inverse of the symplectic transformation that puts $H$ in symplectic normal form.

   Then
   \begin{align}
	   \frac{1}{2} \hat{r}^T H \hat{r} &= \frac{1}{2} \hat{r}^T S_H^T \left( \Lambda \otimes I_2 \right) S_H \hat{r} \\
									   &= \frac{1}{2} \left( S_H \hat{r} \right)^T \left( \Lambda \otimes I_2 \right) S_H \hat{r}.
   \end{align}

   We can then think of $S_H$ as a coordinate transformation and can define a new set of quadrature operators as $ \hat{r}' = S_H \hat{r}$. This is possible since $[ \hat{r}', \hat{r}'^T] = i \Omega$.

   Then the Hamiltonian
   \begin{equation}
	   \frac{1}{2} \hat{r}^T H \hat{r} = \frac{1}{2} \hat{r}'^T (\Lambda \otimes I_2 ) \hat{r}'
   \end{equation}
   is diagonal with respect to this new notation.

   Now consider that
   \begin{align}
	   \frac{1}{2} \hat{r}'^T \left( \Lambda\otimes I_2\right) \hat{r}' &= \frac{1}{2} \sum_{jk} \hat{r}'_{j} \left( \Lambda \otimes I_2 \right)_{jk} \hat{r}'_k \\
																	 &= \frac{1}{2} \sum_j \lambda_j \big[ {\hat{x}}_{j}^{\prime 2} + {\hat{p}}_{j}^{\prime 2} \big].
   \end{align}

\section{Symplectic decomposition of a positive definite matrix}

   Given a positive definite $M$, how can one compute its symplectic matrix $S$ and its symplectic eigenvalues?

   
   To do the task described above, one need only perform the usual eigendecomposition of the matrix $i \Omega M$. Why does this work? By the Williamson Theorem, it follows that
   \begin{equation}
   		M  = S D S^T
   \end{equation}
   for $S$ symplectic, and where $D= \bigoplus_{j=1}^n d_j \left( \begin{matrix} 1&0\\ 0&1 \end{matrix} \right)$ is the diagonal matrix of symplectic eigenvalues of $M$.

   Then consider that
   \begin{align}
   		i \Omega M &= i \Omega S D S^T \\
				   &= i \Omega S (\underline{D} \otimes  I_2) S^T,
   \end{align}
   where $\underline{D} = \begin{bmatrix} d_1 &  &  \\  & \ddots & \\ & & d_n \end{bmatrix}$.

   Since $S\Omega S^T = \Omega$ (due to $S$ being symplectic), we have the following:
   \begin{align}
		&S \Omega S^T \Omega = \Omega^2 = - I \\
	   \implies& S \Omega S^T \Omega S = -S \\
	   \implies& S^{-1} S \Omega S^T \Omega S  = - S^{-1} S = - I \\
	   \implies& \Omega S^T \Omega S = - I \\
	   \implies& \Omega^T \Omega S^T \Omega S = -\Omega^T \\
	   \implies& S^T \Omega S = \Omega \\
	   \implies& \Omega S = S^{-T} \Omega.
   \end{align}
Then we find that
\begin{align}
	i \Omega S \left( \underline{D} \otimes  I_2 \right) S^T &= i S^{-T} \Omega \left( \underline{D} \otimes  I_2 \right) S^T \\
														 &= i S^{-T} \left(  I_n \otimes \Omega_1 \right) \left( \underline{D} \otimes  I_2 \right) S^T \\
														 &= S^{-T} \left( \underline{D} \otimes i \Omega_1 \right) S^T \\
														 &= S^{-T} \left( \underline{D} \otimes -\sigma_Y \right) S^T.
\end{align}
   The last equality follows from the observation that $i \Omega_1 = \left( \begin{matrix} 0&i \\ -i&0 \end{matrix} \right) = -\sigma_Y$, where $\sigma_Y$ is the usual Pauli matrix.

   Also note that $-\sigma_Y = U (-\sigma_Z) U^{\dag}$ with $U = \frac{1}{\sqrt{2}} \left( \begin{matrix} 1&1 \\ i&-i \end{matrix} \right)$. Continuing our previous series of calculations with this observation, we have
   \begin{align}
		i \Omega S \left( \underline{D} \otimes  I_2 \right) S^T
		&= S^{-T} \left( \underline{D} \otimes -\sigma_Y\right) S^T \\
&= S^{-T} \left( \underline{D} \otimes U(-\sigma_Z)U^{\dag} \right) S^T \\
&= S^{-T} \left(  I_n \otimes U \right) \left( \underline{D} \otimes -\sigma_Z \right) \left(  I_n \otimes U^{\dag} \right) S^T.
   \end{align}

   In the above, let us take $B = S^{-T} \left(  I_n \otimes U \right)$ and thus $B^{-1} = \left(  I_n \otimes U^{\dag} \right) S^T$.

   Also, consider that 
   \begin{equation}
	   \underline{D} \otimes -\sigma_Z = \begin{bmatrix} -d_1 \\  &d_1 \\ &&-d_2 \\ &&& d_2 \\ &&&&\ddots \\ &&&&& -d_n \\ &&&&&&d_n \end{bmatrix} 
   \end{equation}

   Now it is apparent that the above calculations made amount to a diagonalization. Thus we have shown that the usual diagonalization of $i \Omega M$ gives eigenvalues $\{ -d_1, d_2, -d_2, d_2, \ldots, -d_n, d_n\}$, which contains the symplectic eigenvalues of $M$.

   Lastly, the eigenvector matrix is $S^{-T} ( I_n \otimes U)$.
   
\section{Relation between Hamiltonian matrix and covariance matrix}


What is the relationship between the Hamiltonian matrix $H$ and the covariance matrix $\sigma$ for a general Gaussian state?

\begin{lemma}
	For a general Gaussian state, the Hamiltonian matrix and covariance matrix are related to each other by
\begin{align}
	\sigma &= \coth\left( \frac{i \Omega H}{2} \right) i \Omega, \\
	H &= 2 \operatorname{arcoth} (i\Omega \sigma) i \Omega,
\end{align}
which can be seen as a generalization of what was found for the diagonal case. 
\end{lemma}

\begin{proof}
Start with positive definite matrix $H$. The symplectic diagonalization is
\begin{equation}
	H = S^T \bigoplus_{j=1}^n \lambda_j \left( \begin{matrix} 1&0 \\ 0&1 \end{matrix} \right) S = S^T D S.
\end{equation}

An earlier argument established that the covariance matrix is
\begin{align}
	\sigma &= S^{-1} \bigoplus_{j=1}^n \coth\left(\frac{\lambda_j}{2}\right) \left( \begin{matrix} 1&0 \\ 0&1 \end{matrix} \right) S^{-T} \\
		   &= S^{-1} \coth\left(\frac{D}{2}\right) S^{-T}
\end{align}

Then consider that, from previous reasoning,
\begin{equation}
	\frac{1}{2} i \Omega H = \frac{1}{2} S^{-1} ( I_n \otimes U)(\underline{D} \otimes -\sigma_Z)( I_n \otimes U^{\dag}) S
\end{equation}
\begin{align}
	\implies &\coth\left(\frac{i\Omega H}{2}\right) \\	
	& = S^{-1} ( I_n \otimes U) \left[\coth\left(\frac{\underline{D} \otimes -\sigma_Z}{2}\right)\right] ( I_n \otimes U^{\dag}) S \\
	& = S^{-1} ( I_n \otimes U) \left[ \coth\left(\frac{\underline{D}}{2}\right) \otimes -\sigma_Z \right] ( I_n \otimes U^{\dag}) S \\
	& = S^{-1} \left[ \coth\left(\frac{\underline{D}}{2}\right) \otimes -\sigma_Y \right]  S \\
	& = S^{-1} \left[ \coth\left(\frac{\underline{D}}{2}\right) \otimes i \Omega_1 \right] S \\
	& = S^{-1} \left[ \coth\left(\frac{\underline{D}}{2}\right) \otimes  I_2 \right] \left[  I_n \otimes i\Omega_1 \right] S \\
	& = S^{-1} \coth\left(\frac{D}{2}\right) i \Omega S \\
	& = S^{-1} \coth\left(\frac{D}{2}\right) S^{-T} i \Omega \\
	& = \sigma i \Omega.
\end{align}

In the above chain of steps, the first equality is due to the functional calculus of matrices. The second equality is due to the $\coth$ function being odd. By simplifying and using properties of $\Omega$, one can simplify towards the end. Summarizing the above, we have shown that
\begin{align}
	 \coth\left(\frac{i \Omega H}{2}\right) = \sigma i \Omega
\end{align}
which implies that
\begin{align}
	\coth\left(\frac{i \Omega H}{2}\right) i \Omega &= \sigma i \Omega (i \Omega) = \sigma.
\end{align}

And thus we are done.

A similar proof can be constructed to yield the reverse result, i.e.
\begin{equation}
	H = 2 \textrm{arcoth} (i \Omega \sigma) i \Omega ,
\end{equation}
concluding the proof.
\end{proof}


%\bibliography{refs}
%\bibliographystyle{alpha}


\end{document}
