In Gaussian Quantum Mechanics, a unitary preserving the Gaussian nature of the state is a called a Gaussian Unitary. In the phase space picture, a Gaussian state is fully characterized by its first and second statistical moments, $\boldsymbol{X}$ and $\sigma$ defined below.
Suppose we have $N$ bosonic modes (or quantum harmonic ocsillators) with the usual commutation relations. Now define the position and momentum operators for the $k^{th}$ mode as, $$\hat{q}_k=\frac{1}{\sqrt{2}}(\hat{a}^\dagger + \hat{a}), \, \hat{p}_k=\frac{-i}{\sqrt{2}}(\hat{a}^\dagger - \hat{a}).$$ Let $\hat{\boldsymbol{X}}=(\hat{q}_1,\hat{p}_1,...,\hat{q}_N,\hat{p}_N)^T,$ then the first statistical moment is the vector of means $\langle \hat{\boldsymbol{X}} \rangle=:\boldsymbol{X}$. The second statistical moment is the covariance matrix $\sigma$ with entries $\sigma_{jk}=\frac{1}{2}\langle \hat{X}_j\hat{X}_k +\hat{X}_k\hat{X}_j\rangle -\langle \hat{X}_j\rangle \langle\hat{X}_k \rangle.$
The most general $\hat{U}_{G}$ Gaussian Unitary acting on a Gaussian state, can equivalently be view in phase space as acting on $\boldsymbol{X}$ and $\sigma$ in the following way: $$ \boldsymbol{X} \mapsto M\boldsymbol{X} + \boldsymbol{\alpha} $$ $$ \sigma \mapsto M\sigma M^T + R$$ for a $2N\times 2N$ real symplectic matrix $M$, a real $2N$ vector $\boldsymbol{\alpha}$ and a $2N \times 2N$ real symmtetric positive semidefinite matrix $R$.
Finally, the covariance matrix can be symplectically diagonalized (it is positive definite and symmetric) as follows. There exists some symplectic matrix $S$ such that $$ \sigma = S \widetilde{\sigma} S^T, \, \text{ where } \widetilde{\sigma}=\bigoplus_{k=1}^N \widetilde{\lambda}_k \mathbb{I}_{2\times 2}. $$
Here $\widetilde{\lambda}_k$ are the symplectic eigenvalues of $\sigma$.
Suppose we have a Gaussian state with covariance matrix $\widetilde{\sigma}$ and the zero vector for its vector of means, we can call its density operator $\hat{\rho}(\vec{0},\widetilde{\sigma})$. Suppose we transform in phase space $\widetilde{\sigma}\mapsto \sigma$, this must correspond to a Gaussian unitary $\hat{U}_S$ such that the new state in the Hilbert space looks like $$\hat{\rho}(\vec{0},\sigma)=\hat{U}_S \hat{\rho}(\vec{0},\widetilde{\sigma})\hat{U}_S^\dagger$$
My question is: given the covariance matrix $\sigma$ and its symplectic eigenvalues $\{\widetilde{\lambda}_k\}$ (or we equivalently have $\widetilde{\sigma}$) and the corresponding sympelctic matrix $S$ taking $\widetilde{\sigma}$ to $\sigma$ through $S\widetilde{\sigma}S^T=\sigma$, how can we find $\hat{U}_S$?
I've looked through the top review articles on Gaussian QM and Quantum Info, they all mention that this correspondence exists but not how to explicitly find the Gaussian unitary from the givens mentioned above.
Any help is appreciated, thanks!
