A general (pure) Gaussian state has the form $\newcommand{\on}[1]{\operatorname{#1}}\newcommand{\ket}[1]{\lvert #1\rangle}\ket{\alpha,\xi}\equiv D(\alpha)S(\xi)\ket{\on{vac}}$, with $\ket{\on{vac}}$ the vacuum state, $D(\alpha)\equiv\exp(\alpha a^\dagger-\alpha^* a)$ the displacement operator, and $S(\xi)\equiv\exp(\frac12 \xi a^{\dagger 2}-\frac12\xi^* a^2)$ the squeezing operator.
I'm interested in the Fock state probabilities of such a state as a function of $\alpha,\xi\in\mathbb C$, that is, the values $$P_k\equiv \langle k|\alpha,\xi \rangle \equiv \frac{1}{\sqrt{k!}}\langle \on{vac}| a^k|\alpha,\xi\rangle, \qquad k\in\mathbb N$$ What are good ways to compute these?
Two "natural" approaches are (1) the direct calculation decomposing the operators as a sum of creations and annihilation operators, and (2) passing through the Wigner representation of Gaussian states. Other "better" routes are likely also viable.
Following the former path, I get, applying BCH and other standard techniques to rearrange the exponential expressions, to the expression: $$ \langle n|\alpha,\xi\rangle = \frac{e^{|\alpha|^2/2}}{\sqrt{n!\cosh(|\xi|)}} \alpha^n \sum_{\ell,k\ge0} (-|\alpha|^2)^\ell (C/\alpha^2)^k \frac{(n+\ell)!}{\ell! k! (\ell+n-2k)!} $$ with $C\equiv \frac12 e^{i\theta_\xi} \tanh(|\xi|)$ and $\xi=|\xi|e^{i\theta_\xi}$. I cannot get this into a useful closed formula. Assuming this is possible at all, knowing the solution obtained from another route (e.g. via the Wigner representation or similar) might shed light into how to handle this.
