How to compute expectation value $\langle e^{iH}\rangle$ for quadratic Hamiltonians? I have a rather basic, but actually non-trivial question:
We consider a bosonic system with creation operators $\hat{a}_i^\dagger$ and annihilation operators $\hat{a}_j$ and vacuum state $|0\rangle$ with $\hat{a}_i|0\rangle=0$. We consider a quadratic Hamiltonian
\begin{align}
\hat{H}=\sum_{ij}(A_{ij}\hat{a}_i\hat{a}_j+A^*_{ij}\hat{a}^\dagger_i\hat{a}^\dagger_j+B_{ij}\hat{a}_i^\dagger\hat{a}_j)\,.
\end{align}
Hermiticity requires $B_{ij}=B^*_{ji}$. How can I compute the expectation value:
\begin{align*}
\langle 0|e^{i\hat{H}}|0\rangle
\end{align*}
In many situations, it is useful to express the problem in terms of Hermitian quadrature operators $\hat{q}_i=\frac{1}{\sqrt{2}}(\hat{a}^\dagger_i+\hat{a}_j)$ and $\hat{p}_i=\frac{1}{\sqrt{2}}(\hat{a}^\dagger_i-\hat{a}_j)$ to define $\hat{\xi}^a=(\hat{q}_1,\dots,\hat{q}_N,\hat{p}_1,\dots,\hat{p}_N,)$, such that $$\hat{H}=\frac{1}{2}h_{ab}\hat{\xi}^a\hat{\xi}^b$$ with implied sums over $a$ and $b$. It is also useful to define the symplectic form $\Omega^{ab}$, such that $[\hat{\xi}^a,\hat{\xi}^b]=i\Omega^{ab}$, which gives rise to the symplectic generator $K^a{}_b=\Omega^{ac}h_{cb}$. We only require that $h_{ab}$ is symmetric and real.
 A: 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$. 
A: 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...
A: Start diagonalizing H. You can follow this:
https://arxiv.org/pdf/0908.0787.pdf
Then use the power series of an exponential. You get the result in the diagonal basis. Make the inverse transformation and it is done.
A: Because I was hoping to get some standard reference, which treats this case, I thought it would be best to keep the question brief.However, as I received a lot of partial suggestions, let me try to explain my thoughts, so far. If you know some standard reference, which treats the general case, it would be greatly appreciated!
For this, let us introduce a basis $\hat{\xi}^a\equiv(\hat{q}_1,\hat{p}_1,\dots,\hat{q}_N,\hat{p}_N)$. With this, we can write
\begin{align}
\hat{H}=\frac{1}{2}h_{ab}\hat{\xi}^a\hat{\xi}^b\,,
\end{align}
where $h_{ab}$ is symmetric and real. It can be constructed from the matrices $A$ and $B$ from the question, but it is useful to write everything in this real basis. Note that $\hat{q}_i=\frac{1}{\sqrt{2}}(\hat{a}_i^\dagger+\hat{a}_i)$ and $\hat{p}_i=\frac{i}{\sqrt{2}}(\hat{a}_i^\dagger-\hat{a}_i)$, such that $\hat{q}_i$ and $\hat{p}_j$ satisfy canonical commutation relations. We can write the commutation relations as $[\hat{\xi}^a,\hat{\xi}^b]=i\Omega^{ab}$, where $\Omega^{ab}$ is a symplectic form.
We can now define the symplectic transformation $M=e^{K}$ with $K^a{}_b=\Omega^{ac}h_{cb}$ which implements the transformation $U=e^{-i\hat{H}}$ as Bogoliubov transformation, namely
\begin{align}
U^\dagger \hat{\xi}^a U=M^a{}_b\hat{\xi}^b\,.
\end{align}
We can further use a polar decomposition (or Cartan decomposition) of $M=Tu$, where $T=\sqrt{MM^\intercal}$, which is well-defined due to $MM^\intercal>0$, and $u=T^{-1}M$. It is further easy to show that $uu^\intercal=1$. We can define $X=\log{T}$ and $Z=\log{u}$. Using the inverse symplectic form $\Omega^{-1}_{ab}$ with $\Omega^{ac}\Omega^{-1}_{cb}=\delta^a{}_b$, we can define two new quadratic Hamiltonians
\begin{align}
\widehat{X}&=\frac{1}{2}\Omega^{-1}_{ac}X^c{}_b\,\hat{\xi}^a\hat{\xi}^b\,,\\
\widehat{Y}&=\frac{1}{2}\Omega^{-1}_{ac}Z^c{}_b\,\hat{\xi}^a\hat{\xi}^b\,.
\end{align}
Using the representation theory of exponential of quadratic Hamiltonians, one can show that we have $e^{-i\hat{H}}=\pm e^{-i\widehat{X}}e^{-i\widehat{Y}}$. Furthermore, one can show that $|0\rangle$ is the ground state of $\widehat{Y}$ and that the expectation value $\langle 0|e^{-i\widehat{X}}|0\rangle$ is real, positive and can be expressed in terms of the spectrum of the matrix $T$, which itself has positive real spectrum (see below). The ground state energy $\widehat{Y}$ is also easily computed from the spectrum of $Y$, which has purely imaginary spectrum, namely the ground state energy is given by $\mathrm{Tr}(\Omega \log{u})/4$.
This gives overall the solution:
\begin{align}
\langle 0| e^{-i\hat{H}}|0\rangle=\pm e^{-i\langle 0|\hat{Y}|0\rangle} \det\left(\frac{\sqrt{2}T^{1/2}}{\sqrt{1+T^2}}\right)=\pm e^{-i\mathrm{Tr}(\Omega \log{u})/4}\det\left(\frac{\sqrt{2}T^{1/2}}{\sqrt{1+T^2}}\right)\,,
\end{align}
where we can bring $T$ into diagonal form $T=\mathrm{diag}(\lambda_1,1/\lambda_1,\dots,\lambda_N,1/\lambda_N)$ with $\lambda_i>0$.
This means my only trouble is the prefactor $\pm$, which is due to the fact that the exponential of quadratic Hamiltonians only form a projective representation of the symplectic group, i.e., $U(M_1)U(M_2)=\pm U(M_1M_2)$.
