How to compute expectation value $\langle e^{iH}\rangle$ for quadratic Hamiltonians?

There is a discussion of the evolution of a wavefunction under your general quadratic hamiltonian in the book by Guillemin and Sternberg "Symplectic methods in Physics." It's in their section on optics. I'm stuck at home because of the lockdown so I do not have my copy to hand. The only thing I have in my notes is the following that I copied from their text. It does discuss the phase issue, but annoyingly with the same sign ambiguity that you have. I only give what I have to give you a sense of what is in G&S, so that you can decide as to whether their book is worth looking up.

Here is their quote: Let $$ M= \left(\matrix{A&B\cr C&D}\right)\in {\rm Sp}(2n, {\mathbb R}) $$ Then the metaplectic projective representation on square integrable functions is given by $$ U(M) \psi(x)=i^\#e^{-i\pi n/4} |\det B|^{-1/2} \frac 1{(2\pi)^{n/2}} \int e^{iW(x,y)}\psi(y) d^ny $$ where $$ W(x,y)=\frac 12 (xDB^{-1} x+ yB^{-1}Ay-2x(B^T)^{-1}y) $$ and $\#$ is even if $\det B>0$ and odd if $\det B<0$.

Using this John Rawnsley: On the universal covering group of the real symplectic group, I found a relatively elegant way to write a partial solution, which makes clear where the sign ambiguity comes from. For this, we use the basis $\hat{\xi}^a=(\hat{q}_1,\dots,\hat{q}_N,\hat{p}_1,\dots,\hat{p}_N)$, such that the ground state covariance matrix $G^{ab}=\langle 0|\hat{\xi}^a\hat{\xi}^b+\hat{\xi}^b\hat{\xi}^a|0\rangle$ is the identity. We then define the associated complex structure and symplectic form \begin{align} G\equiv\begin{pmatrix} 1\!\!1 & 0\\ 0 & -1\!\!1 \end{pmatrix}\,,\quad J\equiv\begin{pmatrix} 0 & 1\!\!1\\ -1\!\!1 & 0 \end{pmatrix}\,,\quad \Omega\equiv\begin{pmatrix} 0 & 1\!\!1\\ -1\!\!1 & 0 \end{pmatrix} \end{align} I use different symbols here, because $J^a{}_b$ is a linear map and $\Omega^{ab}$ is a bilinear form, which only happen to be represented by the same matrix in this specific basis. The quadratic Hamiltonian is fully characterized by its symmetric coefficient matrix $h_{ab}$, such that $\hat{H}=\frac{1}{2}h_{ab}\hat{\xi}^a\hat{\xi}^b$. With this in hand, we define the symplectic algebra generator $K^a{}_b=\Omega^{ab}h_{bc}$ with the symplectic transformation \begin{align} M=e^{t K}=\begin{pmatrix} A & B\\ C & D \end{pmatrix}\,. \end{align} The expectation value is then given by (I derived it with the method of my previous answer based on the definitions of the Rawnsley paper) \begin{align} \langle 0|e^{-i t\hat{H}}|0\rangle=1/\sqrt{\det{\mathcal{C}_M}}\,, \end{align} where $\mathcal{C}_M$ is a complex $N$-by-$N$ matrix defined as follows: We first compute $C_M=\frac{1}{2}(M-JMJ)$, which commutes with $J$, i.e., it satisfies $[C_M,J]=0$. However, $C_M$ is a $2N$-by-$2N$ matrix, but because it commutes with $J$ (representing an imaginary unit with $J^2=-1\!\!1$) we can convert it to complex matrix $N$-by-$N$ matrix in the $(+i)$-eigenspace of $J$. We thus find \begin{align} C_M=\frac{M-JMJ}{2}=\frac{1}{2}\begin{pmatrix} A+D & B-C\\ C-B & A+D \end{pmatrix}\quad\Rightarrow\quad \mathcal{C}_M=\frac{A+D}{2}+i\frac{B-C}{2}\,. \end{align} Now we see exactly where the sign ambiguity comes from: The square root function $1/\sqrt{\det{\mathcal{C}_M}}$ has a branch cut and is thus ambiguous up to a sign. Of course, we can always determine the correct sign by hand: We start at $t=0$ with $\langle 0|e^{-i t\hat{H}}|0\rangle=1$ and then increase $t$ to our desired value, while choosing the square root $1/\sqrt{\det{\mathcal{C}_M}}$ in a continuous way.

Example 1. In the special case, where $[K,J]=0$, we also have $[M,J]=0$, such that $C_M=M=e^{Kt}$ with $A=D$ and $B=-C$. Therefore, we have $\mathcal{C}_M=A+i B$. In this case, $K$ is of the following form and we define $\mathcal{K}$ \begin{align} K=\begin{pmatrix} 0 & k\\ -k & 0 \end{pmatrix}\quad\Rightarrow\quad\mathcal{K}=i k \end{align} where $k=\mathrm{diag}(\epsilon_1,\dots,\epsilon_N)$. We then have $\mathcal{C}_M=A+iB=e^{t\mathcal{K}}=e^{itk}$ as complex $N$-by-$N$ matrix. We can then compute the expectation value as \begin{align} \langle 0|e^{-i t\hat{H}}|0\rangle=1/\sqrt{\det{\mathcal{C}_M}}=1/\sqrt{\det{e^{it k}}}=e^{-\frac{1}{2}\mathrm{tr}\log{ikt}}=e^{-\frac{it}{2}\mathrm{tr}(k)}=e^{-\frac{it}{2}\sum_i\omega_i}\,. \end{align} Here, we got rid of the sign ambiguity by using $\sqrt{\det{e^X}}=e^{\frac{1}{2}\log \det{e^X}}=e^{\frac{1}{2}\mathrm{tr}{\log e^{X}}}$, where we set $\log e^X=X$ without taking branch cuts into account and thereby healing the ambiguity. The result is of course, trivially correct, because for $[J,K]=0$ the state $|0\rangle$ is the ground state of the Hamiltonian with vacuum energy $E_0=\frac{1}{2}\sum_i\omega_i$.

Example 2. The other case that I could treat was the setting where the Hamiltonian $\hat{H}=\sum_i\omega_i(\hat{n}_i+\frac{1}{2})$ has a ground state $|\{0,\dots,0\}\rangle=|0_1\rangle\dots|0_N\rangle$, i.e., it is a tensor product of squeezed 1-mode states, where the initial vacuum $|0\rangle=|0\rangle\dots|0\rangle$ is also a tensor product. Up to a complex phase, we can then relate $|0\rangle$ and $|n_i\rangle$ (arbitrary number excitations of the $n$-th mode of the Hamtiltonian) with a squeezing parameter $r_i$, such that \begin{align} |\langle 2n_i|0\rangle|^2=\frac{(2n)! \tanh^{2n_i}{r_i}}{2^{2n_i}(n_i!)^2\cosh{r_i}}\,. \end{align} This allows us to use a resolution of the identity, sum the series to find \begin{align} \langle 0|e^{-i t\hat{H}}|0\rangle&=\sum_{n_i,n'_j}\langle 0|\{n_1,\dots,n_N\}\rangle\langle\{n_1,\dots,n_N\}|e^{-i \hat{H}}|\{n_1',\dots,n_N'\}\rangle\langle\{n_1',\dots,n_N'\}|0\rangle\\ &=\prod_i\sum_{n_i}|\langle 2n_i|0\rangle|^2e^{-i \omega_i t(2n_i+\frac{1}{2})}\\ &=\prod_i\frac{e^{-\frac{i\omega_it}{2}}\mathrm{sech}(r_i)}{\sqrt{1-e^{-2i\omega_it}\tanh^2{r_i}}}\,. \end{align} However, I haven't really tested if the square root in the denominator has the same sign ambiguity.

I have not been able to come up with a smart way to extract the sign in other cases, except manually taking $t$ to its desired value on requiring continuity of $\langle 0|e^{-i t\hat{H}}|0\rangle$ on the way...

Update: I believe that I found a closed form solution of the problem that does not involve integration. It requires the use of the cocycle function $\eta(M_1,M_2)$ defined in John Rawnsley's paper from above (which can be evaluated with elementary matrix operation). In the end, we find that the correct complex phase is given by $e^{i \psi}$ with $\psi=\frac{1}{2}(\mathrm{tr}(uKu^{-1})+\eta(u,e^K)+\eta(ue^{K},u^{-1}))$, where $u$ is a symplectic transformation that brings $K$ into a standard form $\tilde{K}=uKu^{-1}$ in the same basis as $J$ takes its standard form. More precisely, we transform $K$, such that in its real eigenbasis (for complex eigenvalues we keep real blocks with antisymmetric pieces) the complex structure $J$ takes the same form $\Omega$. I may write this up as a more pedagogical note later on, but for now I hope that my solution is useful if somebody has the same question...