Suppose that I have a single massive quantum mechanical particle in $d$ dimensions ($1\leq d\leq3$), under the action of a well-behaved potential $V(\mathbf r)$, and that I let it settle on the ground state $|\psi_0⟩$ of its hamiltonian, $$ \left[\frac{\mathbf p^2}{2m}+V(\mathbf r)\right]\left|\psi_0\right> = E_0\left|\psi_0\right>. $$ Suppose, because I'm a hater of QM or whatever reason, that I want to model this state of the system as an ensemble of classical trajectories. (More realistically, I might be interested in doing classical trajectory Monte Carlo (CTMC) calculations of some interaction that is hard to model using the full TDSE.) As such, I want to find a classical equivalent to my ground state that I can then use as starting condition for whatever simulation I want to do.

This wavefunction can be examined in phase space in a number of ways, such as using Wigner functions or the Sudarshan $P$ and Husimi $Q$ representations, all of which offer different views into the state and different quasi-classical ways to understand it. My precise question is as follows:

• Is there a way to translate wavefunctions $|\psi⟩$ of a quantum system into probability distributions $\rho(\mathbf r,\mathbf p)$ over classical phase space, in such a way that the eigenstates of a given quantum hamiltonian will be stationary states of the Liouville equation for the corresponding classical system?

That is, I know that after the translation I will get a counterfeit $\rho( \mathbf r, \mathbf p)$ that doesn't actually describe what's going on, but at least I would like it to sit still once I let classical mechanics take over.