Let $\rho$ be the density operator of a Gaussian quantum state on $M$ modes. This implies that its Wigner function can be written as $$ W_{\text{Gaussian}}\left(\boldsymbol{q},\boldsymbol{p}\right)=\frac{1}{\pi^{N}}\frac{1}{\sqrt{\det\left(\sigma\right)}}e^{-\left(\boldsymbol{X}-\boldsymbol{d}\right)^{T}\sigma^{-1}\left(\boldsymbol{X}-\boldsymbol{d}\right)}$$ where $\boldsymbol{X}=(\boldsymbol{q}, \boldsymbol{p})$. Clearly, Gaussian states are fully determined by the covariance matrix $\sigma$ and the displacement vector $\boldsymbol{d}$.
Question: Is there a formula for the $P$ representation of such a state $\rho$ in terms of $\sigma$ and $\boldsymbol{d}$?
I know that the Husimi Q function of a Gaussian state is also a Gaussian, on the other hand, the $P$ function cannot just be a Gaussian since entangled states lead to negativity in the $P$ function.
Edit: Starting a solution: The $P$-function can be obtained from its characteristic function as a Fourier transform (let's start with one mode only for simplicity) $$P(\alpha) = \frac{1}{\pi^2} \int d^2 \lambda \; e^{\lambda^{*}\alpha-\lambda \alpha^{*}} \chi_P(\lambda)$$ and furthermore $\chi_P(\lambda) = \chi_W (\lambda) e^{\frac{1}{2} |\lambda|^2}$ where $\chi_W(\lambda)$ is the characteristic function corresponding to the Wigner function of the state. Hence, $$P(\alpha) = \frac{1}{\pi^2} \int d^2 \lambda \; e^{\frac{1}{2} |\lambda|^2}e^{\lambda^{*}\alpha-\lambda \alpha^{*}} \int d^2\beta \; W(\beta) e^{\lambda\beta^{*}-\lambda^{*} \beta}$$
However, if I try integrating out $\lambda$ now, the integral does not converge if I am not mistaken. So to summarize, I have trouble inverting the known Weierstrass transform $$W(\alpha,\alpha^*)= \frac{2}{\pi} \int P(\beta,\beta^*) e^{-2|\alpha-\beta|^2} \, d^2\beta$$ given on Wikipedia.
