It is known that Gaussian states admit Gaussian Wigner function in phase space, and calculations of their eigenvalues (i.e., diagonalizing the corresponding destiny matrices in Hilbert space) can be done in the Covariant Matrix (CM) level, providing faster way of density matrix diagonalization.
Generally speaking, the reason why we can do that is because Gaussian Wigner function allows normal mode decomposition, and each individual normal mode can be regarded as thermal states of bilinear Hamiltonian, which admits analytic spectrum. There are many excellent books that introduce Gaussian states, the one I like most is https://www.taylorfrancis.com/books/mono/10.1201/9781003250975/quantum-continuous-variables-alessio-serafini.
We know the reason why it is fast for Gaussian states: because we can decompose the density matrix diagonalization in Hilbert space into two steps: the first one is normal mode decomposition in space space, and the second one is Hilbert space diagonalization for each normal mode.
However, Gaussian states are only a subclass of all the quantum states. While we can easily create a manually crafted artificial states that can be normal mode decomposed, it is clear that the above procedure does not apply to all quantum states in phase space. As an example, we know that Wigner function for Fock states involve polynomial terms beyond the phase space Gaussian distribution, which in general cannot be decomposed into multiplications of functions through symplectic transformations on canonical operators.
All the above leads to my current question: is there any ways to do fast eigenvalue calculations beyond the Gaussian states? Hope there is a procedure such that its speed scales with Gaussianity of a state in phase space.
