For boson modes $a^\dagger_i, a_i$, consider the density matrix which is an exponential of quadratic operators in $a^\dagger_i$ and $a_i$: $$\rho = e^{-H_{ij} a^\dagger_i a_j + (K_{ij} a^\dagger_i a^\dagger_j + h.c.)}.$$ Here, $H_{ij}$ and $K_{ij}$ are matrices and Einstein summation convention is assumed.
Question: How can one prove the Wick's theorem $$\langle O X_1 \cdots X_{2n}\rangle = \sum_{\sigma\in S_{2n} \textrm{ contraction}}\langle O X_{\sigma(1)} X_{\sigma(2)} \rangle \cdots \langle O X_{\sigma(2n-1)} X_{\sigma(2n)} \rangle $$ for the above density matrix, where $X_i$ is linear in boson operators and $O$ is an arbitrary operator ordering?
