Can I swap quantum mechanical ground state for some classical trajectory distribution and have it sit still after the swap? 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?


To be fully explicit, I want a map such that produces a classical density with the correct position and momentum distributions, i.e. $\int  \rho(\mathbf r,\mathbf p)\mathrm d\mathbf p = |⟨\mathbf r|\psi⟩|^2$ and $\int  \rho(\mathbf r,\mathbf p)\mathrm d\mathbf r = |⟨\mathbf p|\psi⟩|^2$, and ideally also for all possible quadratures at any angle. The classical $\rho(\mathbf r,\mathbf p)$ should remain stationary under Liouville's equation with a classical hamiltonian $H(\mathbf r,\mathbf p)$ that connects to the quantum hamiltonian via a classical limit or canonical quantization, in the general case, but I'm happy to restrict this to hamiltonians of the form $H(\mathbf r,\mathbf p) = \frac{\mathbf p^2}{2m}+V(\mathbf r)$, in which the correspondence is obvious.
More intuitively, 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.
 A: It seems that OP's sought-for semiclassical realization is possible with the help of the Groenewold-Moyal star product $\star$, i.e. Liouville's eq. becomes
$$ 0~=~\frac{d\rho }{dt}~=~ \frac{1}{i\hbar} [\rho \stackrel{\star}{,}H]+\frac{\partial \rho }{\partial t} .\tag{1}$$ 
In particular, an eigenstate $|\psi\rangle$ of a given quantum Hamiltonian operator $\hat{H}$ (so that the density operator is $\hat{\rho}= |\psi\rangle  \langle \psi|$) will be a stationary state of Liouville's equation (1) for the corresponding semiclassical system, i.e. the Wigner quasiprobability distribution $\rho$ and the Hamiltonian function $H$ will star commute.
We have used the following notation. The $\star$-commutator is defined as
$$ [f\stackrel{\star}{,} g]~:=~f\star g-g\star f.\tag{2}$$ 
The operator-to-function/symbol map $\hat{f}\mapsto f$ is given by the Wigner transform. The inverse map $f\mapsto \hat{f}$ is given by the Weyl transform, cf. e.g. this Phys.SE post. Note that the density operator $\hat{\rho}$ and the Wigner quasiprobability distribution $\rho$ is an example of an operator-symbol/function pair. We stress that that the function/symbol $f$ will in general depend on Planck's constant $\hbar$. See also e.g. my Phys.SE answers here and here. E.g. Heisenberg's operator eq. of motion
$$ \frac{d\hat{f}}{dt}~=~ \frac{1}{i\hbar} [\hat{f},\hat{H}]+\frac{\partial \hat{f}}{\partial t} \tag{3}$$ 
corresponds to
$$ \frac{df}{dt}~=~ \frac{1}{i\hbar} [f\stackrel{\star}{,}H]+\frac{\partial f}{\partial t} \tag{4}$$ 
for functions/symbols $f$.  Liouville's eq. (1) is a special case of eq. (4) with $f=\rho$.
It is possible to generalize the construction to other quasiprobability distributions and operator orderings with their corresponding associative star products. We believe that the use of star products is needed in the generic construction. 
A: Bohm-De-Broglie pilot-wave theory offers a construction of a quantum $\mathcal{O}(\hbar^2)$ correction to the original Hamiltonian along with a classical particle ensemble which does almost exactly what you ask for. The corrected Hamiltonian $H_B$ reads
$$H_B = \frac{\mathbf{p}^2}{2m} + V(\mathbf{x}) - \frac{\hbar^2}{2m} \frac{\Delta \sqrt{\rho}}{\sqrt{\rho}}\,,$$
where $\rho = \int f dp$. The phase-space distribution in $n$ dimensions is then constructed as
$$f(\mathbf{x},\mathbf{p}) = |\psi|^2 \delta^{(n)}(\mathbf{p} - \hbar \Im(\frac{\nabla \psi}{\psi}))\,,$$
where $\psi$ is the original wave-function we want to mimic. It is obvious that at least for the initial condition we will have $\rho=|\psi|^2$ and $$\langle \mathbf{v} \rangle \equiv \int \frac{\mathbf{p}}{m} f dp = \frac{\hbar}{m} \Im(\frac{\nabla \psi}{\psi})\,.$$
What is less obvious is that this ensemble with the Hamiltonian $H_B$ is fully equivalent to the evolution of $\psi$. In other words, the formula for $f$ Liouville-evolved by $H_B$ will be equivalent to the one above in terms of $\psi$ evolved by the Schrödinger equation. This applies not only to ground states but also to general non-stationary states.
One can derive the equivalence by considering the Schrödinger equation, putting $\psi = \sqrt{\rho} \exp(iS/\hbar)$, and interpreting $S$ as the classical Hamilton-Jacobi action. (More details are on the Bohm-De-Broglie wiki page.)
By the way, it seems that something similar to what you are trying to do has been done in quantum-chemistry numerical computations, see the book Applied Bohmian Mechanics: From Nanoscale Systems to Cosmology.
A: If you take the Liouville equation and set $\frac{\partial \rho}{\partial t} = 0$, so that the probability density doesn't depend on the time, you get (in one dimension):
$$\frac{\partial \rho}{\partial r}\dot{r} +\frac{\partial \rho}{\partial p} \dot{p} = 0$$ 
Now, if you use Hamilton equations to know $\dot{r}$ and $\dot{p}$ in terms of the derivatives of the Hamiltonian
$$\frac{\partial \rho}{\partial r}\frac{\partial H}{\partial p} - \frac{\partial \rho}{\partial p}\frac{\partial H}{dr} = 0$$
Which is the Poisson bracket of $\rho$ and $H$. If the probability density is a function of the Hamiltonian only, the Poisson bracket vanishes.
A: I mean, here's how I'd approach it in one dimension...
Consider the coherent states $\alpha~|\alpha\rangle = \hat a |\alpha\rangle$ generated by $\hat x = \lambda (\hat a^\dagger + \hat a)/2,$ and $\hat p = i~\mu(\hat a^\dagger - \hat a)/2.$ They therefore have relatively well defined positions $\langle x \rangle = \lambda \Re(\alpha),$ and momenta $\langle p \rangle = \mu \Im(\alpha);$ in fact they are I believe Gaussians in $(x, p)$ space with only their zero-point fluctuation size. More importantly they resolve the identity with some kernel, $\hat 1 = \int_\mathbb C d^2\alpha~\kappa(\alpha)~ |\alpha\rangle\langle\alpha|.$ Note that even though these are derived for the harmonic oscillator Hamiltonian $\epsilon \hat a^\dagger \hat a,$ which causes them to describe nice circles in phase space just like an actual Harmonic oscillator does, there's no requirement that you actually use them with just that Hamiltonian.
A candidate distribution is therefore to represent $|\psi_0\rangle$ as a function $\mathbb C \to \mathbb C$ as $\psi_0(\alpha) = \langle\alpha|\psi_0\rangle,$ from which we might recover some density $\rho(\alpha) = \kappa(\alpha)~\psi_0^*(\alpha) ~ \psi_0(\alpha),$ losing the quantum phase but necessarily preserving the $\int~d^2\alpha ~\rho(\alpha) = 1.$ Then there is a nice interpretation of $\rho(x, p) = \rho(x/\lambda + i p/\mu)/(\lambda \mu),$ essentially running the above expectation values "in reverse".
The integration $\int dp~\rho(x, p) = \iint dx' ~dp~ \delta(x-x')~\rho(x',p)$ then becomes:$$\int_{\mathbb C} d^2\alpha~\kappa(\alpha)\langle \psi_0|\delta(x - \Re\alpha)|\alpha\rangle\langle\alpha|\psi_0\rangle.$$I'd need to think a little about whether we can replace that $\delta(x - \Re \alpha)$ with the observable $|x\rangle\langle x|,$ but assuming that this is true, the resolution of the identity does all of the rest of the work for us and we'd simply get $\langle\psi_0|x\rangle\langle x|\psi_0\rangle$ as requested, and there's no obvious reason for this to be specific to any one quadrature than to hold up for for all quadratures. So it's a very "obvious" way to construct something which is "mostly" right.
The last thing to prove is that the resulting $\rho(\alpha)$ also is stationary under the Liouville's equation, but given that it comes from a $|\psi_0\rangle$ which is stationary under the Hamiltonian, it seems likely, as a special case of Ehrenfest's theorem... so it's a simpleminded construction but I wouldn't rule it out.
A: I am late in the discussion, and I've missed out on subtleties; but the way I  might approach the problem would be through simple paradigms: well, the oscillator, is as classical as quantum systems go. With proper normalizations, the trajectory of every system for it is rigid rotation in phase space, both for any Wigner function (Groenewold, 1946) and its classical limit Liouville density: $f(x,p;t)=f(x\cos t-p\sin t,p\cos t+x\sin t;0)$ ,

A stationary quantum system is an axially symmetric Wigner function, and same for its classical Liouville density: What happens is that a stationary configuration is but the integral over all phases (call them starting times!) of any configuration. Since everything rotates in lockstep, the configuration appears (is) stationary, both quantum and classical mechanically. 
For instance, this is the 2nd excited state of the oscillator, at scales $O(\hbar)$, 

but there is nothing wrong with contemplating a huge (macroscopic) 
collective system in phase space ... of freshman lab size, which is axially symmetric, a thick swarm of noninteracting charged particles oscillating in a field in lockstep having started their oscillations uniformly at all moments of the oscillation cycle. (You did't say "realistic" now...)
The initial distribution may be picked to be positive semidefinite and axially symmetric: it need not be a pure star-eigenstate: It will rotate rigidly nonetheless and simulate a Liouville density. 
But, in any case, this is the easiest way to illustrate the concept: Stationarity, at least for periodic motion, may well amount to a perfectly phase-smeared ensemble.
The oscillator is, however, in fact, special, to the extent that  the Wigner density evolute is the very same function of the evolved x and p and an axially symmetric configuration is stationary. Virtually unique, except for systems mappable to it. This is just a conceptual existence demonstration, not a method. I am not sure G Braunss 2009 is useful, but there it is...
