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.

  • $\begingroup$ I'm sure I am way off the mark here but I designed a program to simulate 2D snap shots of slit fringe patterns based on classical trajectories. You can find some of my simulators at billalsept.com $\endgroup$ Aug 17, 2016 at 15:53
  • 1
    $\begingroup$ Sorry, Bill, but I don't see any relevance at all of your personal theories to this question, and I don't think this is the appropriate place to advertise them. $\endgroup$ Aug 17, 2016 at 16:01
  • $\begingroup$ It's not a theory it's a simulator and I was just trying to help. I said I might be off the mark. It sounded like you wanted to derive a particle or classical trajectory solution. $\endgroup$ Aug 17, 2016 at 16:20
  • $\begingroup$ A couple of ideas: 1) Based on arxiv.org/abs/0810.2394, see Eq.(10) therein, start as usual with a polar decomposition $\psi=\sqrt{\rho}e^{iS/\hbar}$, identify ${\bf p}\sim\nabla S$, and define a "phase-space probability distribution" as $\rho({\bf x}, t)\delta({\bf p}-\nabla S)$. Then set up the Lagrangian, Hamiltonian etc as discussed at length in paper. 2) Sec.IV.2 in arxiv.org/abs/0712.1984 gives another way to equiv. Hamilton-Jacobi setup, maybe it's worth a look. Once H-J eqs. are in place, translating to/from Schroedinger should be straightforward. $\endgroup$
    – udrv
    Aug 20, 2016 at 21:12
  • 1
    $\begingroup$ Comment to the post (v1): By the word probability distribution, do you mean quasi-probability distribution? $\endgroup$
    – Qmechanic
    Aug 31, 2016 at 18:29

5 Answers 5


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.


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.

  • $\begingroup$ Hang on, I'm not sure I get this. If I understand it correctly, you're saying that the Wigner distribution already follows a Liouville-like equation, except that it needs to be corrected with higher powers of $\hbar$? $\endgroup$ Aug 31, 2016 at 21:57
  • 1
    $\begingroup$ $\uparrow$ Yes. $\endgroup$
    – Qmechanic
    Aug 31, 2016 at 22:33
  • $\begingroup$ Hmmm, ok. To bring this back to the Classical Trajectory Monte Carlo (CTMC) spirit, then, can you comment on whether the modified Liouville equation admits solution by characteristics, and what the equivalent of Hamilton's equations would be for those trajectories? $\endgroup$ Aug 31, 2016 at 22:43
  • $\begingroup$ The method of quantum characteristics is discussed on this Wikipedia page. $\endgroup$
    – Qmechanic
    Sep 1, 2016 at 14:10

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.

  • $\begingroup$ That's definitely interesting. Two comments: (i) it's not clear how restrictive the Ansatz is - my guess is 'quite a bit'. More importantly, (ii) the quantum state has sort of vanished, though it can be reinstated through its eigenvalue as $\rho(\mathbf r,\mathbf p) = \delta(H(\mathbf r,\mathbf p) - E_n)$. What guarantees are there, then, that the resulting marginals on position and momentum will match the quantum eigenstates? If this fails, then it's interesting but not as useful for CTMC applications. $\endgroup$ Aug 30, 2016 at 0:34
  • 1
    $\begingroup$ That said, there's a quicker way to get there - simply noting that the stationary Liouville equation is of the form $\{\rho, H\}=0$, with $\{·,·\}$ the Poisson bracket, and that any $\rho=f(H)$ will trivially satisfy that. There probably is, in fact, a converse showing that form from the equation - but whether it's useful or not for connecting with the quantum eigenstates is not that trivial. $\endgroup$ Aug 30, 2016 at 0:37
  • $\begingroup$ As a idea, if you take the classical probability distribution $\rho = 1/\dot{r}(r)$, with $\dot{r} = \sqrt{2m(E -V(r))}$, it looks like the quantum distribution if the quantum numbers are big. (See for example this image ). This probability density might be stationary or near stationary in classical mechanics and the quantum eigenstates' probability density converges to it. $\endgroup$
    – lytex
    Sep 1, 2016 at 9:58

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.

  • $\begingroup$ Sure, it looks reasonable, but it's not at all obvious to me that that density will be stationary. You could say the same about the Husimi Q and the Sudarshan P representations (they're both baked from $|\psi⟩$ combined with $|\alpha⟩$), but if the same procedure works perfectly for all three then there's definitely something nontrivial going on there. $\endgroup$ Aug 31, 2016 at 22:09
  • $\begingroup$ Hm. Peeking inside and ignoring $\lambda,\mu$ for a second gives $|x+i p\rangle=\exp(-(x^2+p^2)/2)\sum(x+ip)^n/\sqrt{n!}~|n\rangle$ where the $|n\rangle$ are SHO eigenfunctions; this seems to prove $\partial_x|x + i p\rangle =-x |x+ip\rangle + a^\dagger a/(x + i p) |x + i p\rangle,$ and the first $x$ can be of course replaced with $(a + a^\dagger)/2$ or whatever it is. Then this gets multiplied with $\partial H/\partial p,$ but if we did the same trick to $\hat H$ to get $H(x,p) = \langle x+ip|\hat H|x + i p\rangle$ we'd have a further set of terms. I'll definitely think about this more. $\endgroup$
    – CR Drost
    Sep 1, 2016 at 18:33

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)$ , enter image description here

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)$, enter image description here

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...

  • $\begingroup$ This is somewhat problematic if the Wigner function is negative, but more importantly I suspect it only works because of the isoperiodicity of the harmonic oscillator. There you can make a stationary state by starting off with some density and then averaging over one period, but if you have an anharmonic system, over what time interval do you average? On the other hand, if this can be made to work e.g. for a 1D anharmonic oscillator, or for a hydrogenic problem (with definite $l$, say) then it would be much stronger. $\endgroup$ Mar 24, 2017 at 9:02
  • $\begingroup$ As you say, the negative values qualm is dismissible: You may handpick your positive-definite initial symmetric configuration, a mixed state, not a stargenstate, to be axially symmetric, and it stays that way. Indeed, I'm only illustrating/intuiting a principle for the oscillator, which is almost classical, and any other system mappable to it, by hook or crook.... $\endgroup$ Mar 24, 2017 at 10:40
  • $\begingroup$ Sure, but the main thrust of the question is looking at some arbitrary potential, which may be very far from harmonic, and making do with whatever's possible, rather than try to restrict attention to whatever can be shoehorned into the one case that works perfectly. $\endgroup$ Mar 24, 2017 at 10:47
  • $\begingroup$ Oh, you want a method... not an existence demonstration guide. Never seen one. But, isn't this just why we have QM? $\endgroup$ Mar 24, 2017 at 13:12
  • $\begingroup$ Yes, that is why we have QM ;-), but sometimes the TDSE is very expensive (such as the kinds of TDSE solutions by A Scrinzi, for flavour), so people often recur to Classical Trajectory Monte Carlo methods (example) to get a bit more insight into dynamical processes without breaking the bank. This question asked whether you can get a reasonable starting point for the ground state for that kind of simulation (where your observable is e.g. ionization with rescattering). $\endgroup$ Mar 24, 2017 at 13:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.