Quantum flux and reverse engineering of quantum wave functions

The probability flux, or probability current, is introduced in quantum mechanics textbooks as a well-defined operator that depends on position and its expectation value can be calculated exactly at every point in space without uncertainty bounds, but its connection to experiment is often left to the reader's imagination. The flux operator, whose expectation over the wave function gives the traditional flux ${\mathbf { j}}\left ({\mathbf { r}},{\mathbf { p}}\right )$, is defined as

$$\begin{equation} \hat{\mathbf{j}}_{\mathbf{r}}=\frac{1}{2m}\left(\left|\mathbf{r}\vphantom{\mathbf{r}}\right\rangle \left\langle \vphantom{\mathbf{r}}\mathbf{r}\right|\hat{{\bf p}}+\hat{{\bf p}}\left|{\bf r}\vphantom{{\bf r}}\right\rangle \left\langle \vphantom{{\bf r}}{\bf r}\right|\right). \end{equation} \tag{ 1 }$$

Here r and p indicate the position and momentum of a particle while m is the mass. The concept of "flux at a point" seems paradoxical because we say something about momentum while also knowing position precisely. This raises the question: Can the flux even be measured?

On the other hand, probability flux vanishes on stationary states for systems with time-reversal symmetry. This is a shame, since strong semiclassical connections between trajectory flow and quantum eigenstates lie completely hidden in the universal value of 0 for the flux. This letter addresses the problem of the underlying dynamics encoded in stationary states by extending the definition of the flux to coherent-state projections and using a "processed Husimi" representation which makes decisions using many Husimi projections at each location and spits out classical rays. Not only does our approach resolve the measurement problem for the flux, we also use it as a technique for extracting semiclassical paths from a quantum wave function even when the flux is zero.

Several discussions connecting the flux to experimental measurement [1–4] exist in the literature; we begin our argument instead by identifying the eigenstates of the flux operator [5–7], and present a physical interpretation. The flux eigenstates had previously been studied in the context of thermal rates [5, 7] and of wave packets [6]. By providing a clear derivation and a novel application of the flux eigenstates, we provide a perspective of how these techniques can fit into the broader context of wave function analysis.

We begin by replacing the Dirac basis implicit in eq. (1) with the Gaussian basis defined as

$$\begin{equation} \left\langle \left.{\bf r}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\right\rangle =N_{\sigma}^{d/2}e^{-\left({\bf r}-{\bf r}_{0}\right)^{2}/4\sigma^{2}}, \end{equation} \tag{ 2 }$$

where d is the number of dimensions in the system and $N_{\sigma }=\left (\sigma \sqrt {2\pi }\right )^{-1}$ is a normalization constant. There is certainly no loss of generality here; the $\sigma \rightarrow 0$ Dirac delta basis limit is always just a step away. Instead, there is new capability introduced, as we show below. The Gaussian basis over all centers is overcomplete, so nothing is left hidden from the analysis we develop now. Below we are led to a harmonic-oscillator basis at each location. In the Gaussian basis, the modified flux operator is

$$\begin{eqnarray*} &&\hat{{\bf j}}_{{\bf r}_{0},\sigma}=\frac{1}{2m}\left(\left|{\bf r}_{0},\sigma\vphantom{{\bf r}_{0},\sigma}\right\rangle \left\langle \vphantom{{\bf r}_{0},\sigma}{\bf r}_{0},\sigma\right|\hat{{\bf p}}+\hat{{\bf p}}\left|{\bf r}_{0},\sigma\vphantom{{\bf r}_{0},\sigma}\right\rangle \left\langle \vphantom{{\bf r}_{0},\sigma}{\bf r}_{0},\sigma\right|\right). \end{eqnarray*}$$

The eigenstates, projected onto each orthogonal spatial dimension i, are obtained using the eigenvalue equation

$$\begin{equation} \hat{j}_{{\bf r}_{0},\sigma,i}\left|\lambda_{\sigma,i}\right\rangle =\lambda_{\sigma,i}\left|\lambda_{\sigma,i}\right\rangle , \end{equation} \tag{ 3 }$$

which has a solution of the form

$$\begin{eqnarray*} &&\left|\lambda_{\sigma,i}\right\rangle =\left|{\bf r}_{0},\sigma\right\rangle +a\hat{p}_{i}\left|{\bf r}_{0},\sigma\right\rangle . \end{eqnarray*}$$

Using the following two equations,

$$\begin{eqnarray*} &&\left\langle {\bf r}\vphantom{\hat{{\bf p}}}\vphantom{{\bf r}_{0},\sigma}\left|\hat{{\bf p}}\vphantom{{\bf r}}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\vphantom{{\bf r}}\vphantom{\hat{{\bf p}}}\right\rangle =i\hbar\sigma^{-2}\left({\bf r}-{\bf r}_{0}\right)e^{-\left({\bf r}-{\bf r}_{0}\right)^{2}/4\sigma^{2}}\end{eqnarray*}{}$$

and $\left \langle {\mathbf { r}}_{0},\sigma \vphantom {\hat {{\mathbf { p}}}}\vphantom {{\mathbf { r}}_{0,}\sigma }\left |\hat {{\mathbf { p}}}\vphantom {{\mathbf { r}}_{0},\sigma }\vphantom {{\mathbf { r}}_{0,}\sigma }\right |{\mathbf { r}}_{0,}\sigma \vphantom {{\mathbf { r}}_{0},\sigma }\vphantom {\hat {{\mathbf { p}}}}\right \rangle =0,$ we can write

$$\begin{equation} \hat{j}_{{\bf r}_{0},\sigma,i}\left|\lambda_{\sigma,i}\right\rangle =\frac{1}{2m}\left(a\left\langle \hat{p}_{i}^{2}\right\rangle _{\sigma}\left|{\bf r}_{0},\sigma\right\rangle +\hat{p}_{i}\left|{\bf r}_{0},\sigma\right\rangle \right). \end{equation}{} \tag{ 4 }$$

Finding the conditions on λ_σ,i that allow eq. (4) to be written in the form of eq. (3), we obtain

$$\begin{eqnarray*} &&\lambda_{\sigma,i}=\frac{a}{2m}\left\langle \hat{p}_{i}^{2}\right\rangle _{\sigma};\quad\lambda_{\sigma,i}=\frac{1}{2ma}. \end{eqnarray*}$$

Since $\left \langle \hat {p}_{i}^{2}\right \rangle _{\sigma }=\frac {\hbar ^{2}}{4\sigma ^{2}}$, we find the value of $a=\pm \frac {2\sigma }{\hbar }$ from which we obtain the two eigenvalues

$$\begin{equation} \lambda_{\sigma,i,\pm}=\pm\frac{\hbar}{4m\sigma}. \end{equation} \tag{ 5 }$$

The eigenstates take the form

$$\begin{equation} \left\langle \left.{\bf r}\vphantom{\lambda_{\sigma,i,\pm}}\right|\lambda_{\sigma,i,\pm}\right\rangle =\left\langle \left.{\bf r}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\right\rangle \pm\frac{i}{\sigma}{\bf e}_{i}\cdot\left({\bf r}-{\bf r}_{0}\right)\left\langle \left.{\bf r}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\right\rangle , \end{equation} \tag{ 6 }$$

where e_i is the unit vector along spatial direction i. Equation (6) is a linear combination of two functions: the Gaussian (eq. (2)) and its derivative. Projection of a wave function onto the first term can be interpreted as measuring its probability amplitude at point r₀, and projection onto second term as measuring its derivative along the i-th direction at the point r₀.

To determine the expectation value of the flux operator, we begin by labeling the excited states of the harmonic oscillator at position r₀ oriented along the i-th direction,

$$\begin{eqnarray*} \left\langle \left.{\bf r}\vphantom{0}\right|0\right\rangle & = & \left\langle \left.{\bf r}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\right\rangle\!,\\ \left\langle \left.{\bf r}\vphantom{1}\right|1\right\rangle & = & \frac{{\bf e}_{i}\cdot({\bf r}-{\bf r}_{0})}{\sigma}\left\langle \left.{\bf r}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\right\rangle\!,\\ \left\langle \left.{\bf r}\vphantom{2}\right|2\right\rangle & = & \sqrt{\frac{1}{2}}\left(\frac{\left({\bf e}_{i}\cdot({\bf r}-{\bf r}_{0})\right)^{2}}{\sigma^{2}}-1\right)\left\langle \left.{\bf r}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\right\rangle\!;\;\ \mbox{etc.} \end{eqnarray*}$$

These states form a complete set in which the flux operator can be explicitly expressed as the Hermitian matrix

$$\begin{eqnarray*} \hat{j}_{{\bf r}_{0},\sigma,i}=\left(\begin{array}{@{}c@{\qquad}c@{\qquad}c@{\qquad}c@{\qquad}c@{}} 0 & +i\lambda & 0 & \cdots & 0\\ -i\lambda & 0 & 0 & \cdots & 0\\ 0 & 0 & 0 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & 0 \end{array}\right), \end{eqnarray*}$$

where $\lambda =\lambda _{\sigma ,i,+}=\frac {\hbar }{4m\sigma }$. There are additional sets of harmonic-oscillator states centered at points other than r₀ also with zero components in the flux matrix.

The complete set of eigenstates $\left |\lambda _{1}\right \rangle ,\left |\lambda _{2}\right \rangle ,\left |\lambda _{3}\right \rangle ,\dots$ of the flux operator expressed in terms of excited states of the harmonic oscillator are

$$\begin{eqnarray*} &&\left|\lambda_{1}\right\rangle ,\left|\lambda_{2}\right\rangle ,\left|\lambda_{3}\right\rangle ,\dots=\left(\begin{array}{c} +1\\ -i\\ 0\\ \vdots \end{array}\right),\left(\begin{array}{c} +1\\ +i\\ 0\\ \vdots \end{array}\right),\left(\begin{array}{c} 0\\ 0\\ 1\\ \vdots \end{array}\right),\cdots \end{eqnarray*}$$

with eigenvalues −λ, λ, and 0. Measurement by the flux operator collapses the wave function onto one of these eigenstates, the infinite majority of which are in the degenerate zero-eigenvalue subspace spanning all excited states of the harmonic oscillator above |1〉. Only the first two eigenstates yield non-zero flux values. In the $\sigma \rightarrow 0$ limit, the eigenvalues of these two states tend towards positive and negative infinity.

When expanding the flux expectation value, we can use the complete eigenbasis to show that

$$\begin{eqnarray} \langle \psi\vphantom{\hat{j}_{{\bf r}_{0},\sigma,i}}\vphantom{\psi}|\hat{j}_{{\bf r}_{0},\sigma,i}\vphantom{\psi}\vphantom{\psi}|\psi\vphantom{\psi}\vphantom{\hat{j}_{{\bf r}_{0},\sigma,i}}\rangle & = & \left\langle \psi\vphantom{\hat{j}_{{\bf r}_{0},\sigma,i}\sum_{i=1}^{\infty}\left|\vphantom{\sum_{i=1}^{\infty}}\lambda_{i}\vphantom{\lambda_{i}}\right\rangle \left\langle \vphantom{\vphantom{\sum_{i=1}^{\infty}}\lambda_{i}}\lambda_{i}\right|}\vphantom{\psi}\left|\hat{j}_{{\bf r}_{0},\sigma,i}\sum_{i=1}^{\infty}\left|\vphantom{\sum_{i=1}^{\infty}}\lambda_{i}\vphantom{\lambda_{i}}\right\rangle\right.\right. \langle \lambda_{i}|\vphantom{\psi}\vphantom{\psi}|\psi\rangle \nonumber & = & \lambda\left|\left\langle \left.\psi\vphantom{\lambda_{1}}\right|\lambda_{1}\right\rangle \right|^{2}-\lambda\left|\left\langle \left.\psi\vphantom{\lambda_{2}}\right|\lambda_{2}\right\rangle \right|^{2}. \end{eqnarray} \tag{ 7 }$$

From eq. (6), it can be shown that some contributions from $\left |{\left \langle {\psi }|0\right \rangle }\right |^{2}$ and $\left |{\left \langle {\psi }|1\right \rangle }\right |^{2}$ cancel themselves due to the opposite sign of the eigenvalues, and only the cross-term ${\left \langle {\psi }|0\right \rangle }^{\ast }{\left \langle {\psi }|1\right \rangle }-{\left \langle {\psi }|0\right \rangle }{\left \langle {\psi }|1\right \rangle }^{\ast }$ remains. This form is directly related to the following expression for the flux at point r₀:

$$\begin{eqnarray*} &&{\bf j}_{{\bf r}_{0}}\left(\Psi({\bf r})\right)=\frac{\hbar}{2mi}\left(\Psi^{*}({\bf r}_{0})\mathbf{\nabla}\Psi({\bf r}_{0})-\Psi({\bf r}_{0})\mathbf{\nabla}\Psi^{*}({\bf r}_{0})\right). \end{eqnarray*}$$

Equation (7) can be rewritten as

$$\begin{eqnarray} \left\langle \psi\vphantom{\hat{j}_{{\bf r}_{0},\sigma,i}}\vphantom{\psi}\left|\hat{j}_{{\bf r}_{0},\sigma,i}\vphantom{\psi}\vphantom{\psi}\right|\psi\vphantom{\psi}\vphantom{\hat{j}_{{\bf r}_{0},\sigma,i}}\right\rangle & = & \frac{i\hbar}{4m\sigma^{2}}[\left\langle \psi\vphantom{{\bf e}_{i}\cdot\left({\bf r}-{\bf r}_{0}\right)}\vphantom{{\bf r}_{0},\sigma}\left|{\bf e}_{i}\cdot\left({\bf r}-{\bf r}_{0}\right)\vphantom{\psi}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\vphantom{\psi}\vphantom{{\bf e}_{i}\cdot\left({\bf r}-{\bf r}_{0}\right)}\right\rangle \left\langle \left.\psi\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\right\rangle ^{\ast}\nonumber \\ & & -\left\langle \psi\vphantom{{\bf e}_{i}\cdot\left({\bf r}-{\bf r}_{0}\right)}\vphantom{{\bf r}_{0},\sigma}\left|{\bf e}_{i}\cdot\left({\bf r}-{\bf r}_{0}\right)\vphantom{\psi}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\vphantom{\psi}\vphantom{{\bf e}_{i}\cdot\left({\bf r}-{\bf r}_{0}\right)}\right\rangle ^{\ast}\left\langle \left.\psi\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\right\rangle ]. \end{eqnarray} \tag{ 8 }$$

The traditional flux operator corresponds to the limit $\sigma \rightarrow 0$, at which point the two terms in eq. (6) become the delta function and its derivative, while the flux values of the first two eigenstates become

$$\begin{eqnarray*} &&\lim_{\sigma\rightarrow0^{+}}\lambda_{\sigma,i,\pm}=\pm\infty. \end{eqnarray*}$$

In addition, there are an infinite number of other eigenstates with zero eigenvalues. The infinite value for the two non-zero eigenstates is corroborated by the observation of Park and Light [5]; however, in their work they reach this conclusion in the limit of an infinitely large basis [7].

A single application of the flux operator at a particular point in space almost always results in zero, but occasionally in an immensely positive or negative value. It is thus necessary to perform the averaging over an infinite number of measurements to obtain an expression equivalent to the textbook flux.

The prefactors before the Gaussian states in eq. (6) are related to the Taylor expansion,

$$\begin{eqnarray*} &&e^{\pm\frac{i}{\sigma}{\bf e}_{i}\cdot({\bf r}-{\bf r}_{0})}\approx1\pm\frac{i}{\sigma}{\bf e}_{i}\cdot({\bf r}-{\bf r}_{0}).\end{eqnarray*}$$

This suggests a connection between the flux eigenstates and the coherent state, defined as

$$\begin{equation} \left\langle \left.{\bf r}\vphantom{{\bf r}_{0},{\bf k}_{0},\sigma}\right|{\bf r}_{0},{\bf k}_{0},\sigma\right\rangle =N_{\sigma}^{d/2}e^{-\left({\bf r}-{\bf r}_{0}\right)^{2}/4\sigma^{2}+i{\bf k}_{0}\cdot{\bf r}}, \end{equation} \tag{ 9 }$$

which is a Gaussian envelope over the plane wave e^ik₀·r. Observing that the phase e^ik₀·r₀ is arbitrary, we expand in the limit of k₀σ ≪ 1 to produce

$$\begin{eqnarray} \left\langle \left.{\bf r}\vphantom{{\bf r}_{0},{\bf k}_{0},\sigma}\right|{\bf r}_{0},{\bf k}_{0},\sigma\right\rangle & \approx & N_{\sigma}^{d/2}e^{-\left({\bf r}-{\bf r}_{0}\right)^{2}/4\sigma^{2}}\left(1+i{\bf k}_{0}\cdot\left({\bf r}-{\bf r}_{0}\right)\right)\nonumber \\ & \approx & \left\langle \left.{\bf r}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\right\rangle +i{\bf k}_{0}\cdot\left({\bf r}-{\bf r}_{0}\right)\left\langle \left.{\bf r}\vphantom{{\bf r}_{0},\sigma}\right|{\bf r}_{0},\sigma\right\rangle . \end{eqnarray} \tag{ 10 }$$

Note that the dispersion relation for the free-particle continuum is a circle with radius $k_{0}=\frac {\sqrt {2mE}}{\hbar }$, which does not depend on the orientation of k₀. The second terms in eq. (10) and eq. (6) are proportional to each other when k₀ points along the i-th direction. This similarity allows us to relate the flux expectation value from eqs. (7) and (8) to coherent-state projections as

$$\begin{eqnarray} \lim_{\sigma k_{0}\rightarrow0}\left\langle \psi\vphantom{\hat{j}_{{\bf r}_{0},\sigma,i}}\vphantom{\psi}\left|\hat{j}_{{\bf r}_{0},\sigma,i}\vphantom{\psi}\vphantom{\psi}\right|\psi\vphantom{\psi}\vphantom{\hat{j}_{{\bf r}_{0},\sigma,i}}\right\rangle & = & \frac{\hbar k_{0}}{4m\sigma^{2}}\left[\left|\left\langle \left.\psi\vphantom{{\bf r}_{0},k_{0}{\bf e}_{i},\sigma}\right|{\bf r}_{0},k_{0}{\bf e}_{i},\sigma\right\rangle \right|^{2}\right.\nonumber \\ &&-\left. \left|\left\langle \left.\psi\vphantom{{\bf r}_{0},-k_{0}{\bf e}_{i},\sigma}\right|{\bf r}_{0},-k_{0}{\bf e}_{i},\sigma\right\rangle \right|^{2}\right], \end{eqnarray} \tag{ 11 }$$

where the traditional flux vector is constructed from the components in each direction.

Unlike the Gaussian states appearing in the flux eigenstates in eq. (6), the meaning of coherent states is straightforward: they describe a semiclassical particle minimizing the product of position and momentum uncertainties. Using the uncertainty relation Δx∝1/Δk∝σ in the limit $\sigma \rightarrow 0$, results in coherent-state measurements with infinite uncertainty in k-space and zero uncertainty in real space. This is the limit where the traditional flux operates.

The representation of quantum-mechanical systems in terms of phase-space distribution functions, such as Wigner's quasi-probability distribution function [8, 9] and the Husimi distribution function [10], allows expressing quantum-mechanical results into a form that resembles classical mechanics [11–13]. The standard Husimi projection technique consists of an idealized measurement of position and momentum consistent with the uncertainty principle. The Husimi distribution consists of the expectation value of the projector onto a coherent state and is a phase-space probability density.

Based on the coherent-state projection of the flux operator, we can now develop a deeper intuition. The processed Husimi technique uses coherent states to produce maps of the flux that display information about the local phase-space distribution of the wave function in real space. The current flow maps adhere to the uncertainty principle constraints. Applying the processed Husimi representation to a quantum wave function can reverse engineer or deconstruct it, yielding the underlying classical rays even when those rays cross in several directions, as we demonstrate below.

The properties of coherent states make them a suitable basis for expanding the flux operator to a measurable quantity. It is defined as a measurement of a wave function ψ(r) by a coherent state, or "test wave packet", written as

$$\begin{eqnarray} \textrm{Hu}\left({\bf r}_{0},{\bf k}_{0},\sigma;\psi({\bf r})\right) & = & \left|\left\langle \left.\psi\vphantom{{\bf r}_{0},{\bf k}_{0},\sigma}\right|{\bf r}_{0},{\bf k}_{0},\sigma\right\rangle \right|^{2}. \end{eqnarray} \tag{ 12 }$$

Weighting each of these measurements by the wave vector produces a Husumi vector.

In fig. 1, we present sets of thirty-two equally spaced Husimi vectors for two wave functions. Setting $\sigma \rightarrow 0$ results in coherent-state measurements with infinite uncertainty in k-space and zero uncertainty in real space, the limit where the traditional flux operates. For small σ, due to the large momentum uncertainty, coherent-state projections merely reproduce the probability amplitude $\left |\psi ({\mathbf { r}})\right |^{2}$ in all directions of k₀. The flux emerges as a small residual which can be retrieved by summing each coherent-state projection weighted by k₀, such that

$$\begin{equation} \left\langle \psi\vphantom{\hat{{\bf j}}_{{\bf r}_{0},\sigma}}\vphantom{\psi}\left|\hat{{\bf j}}_{{\bf r}_{0},\sigma}\vphantom{\psi}\vphantom{\psi}\right|\psi\vphantom{\psi}\vphantom{\hat{{\bf j}}_{{\bf r}_{0},\sigma}}\right\rangle \approx\int{\bf k}_{0}\left|\left\langle \left.\psi\vphantom{{\bf r}_{0},{\bf k}_{0},\sigma}\right|{\bf r}_{0},{\bf k}_{0},\sigma\right\rangle \right|^{2}{\rm d}^{d}k_{0}. \end{equation} \tag{ 13 }$$

Note that in the limit $\sigma \rightarrow 0$, the contributing points in the k-space integral reduce to just the orthogonal directions. The absence of flux in time-reversal symmetric states can be interpreted as the mutual cancellation of coherent-state projections along each direction in k-space.

Zoom In Zoom Out Reset image size

For larger σ, reduced momentum uncertainty allows for substantial variation in the coherent-state projections as the direction of k₀ varies. Away from the $\sigma \rightarrow 0$ regime, we can use coherent states (Husimi projections) to produce a map of the local phase space of a wave function. This can be seen in fig. 1: as uncertainty is reduced, the uniform sunburst (d) contracts into lobe(s) (c) and finally to unambiguous vector(s) (a and b). By taking snapshots of the phase space at many points across a system for larger σ, we can process the result to produce a semiclassical map showing the dominant classical paths contributing to a given wave function. Thus, the term processed Husimi for these visualizations. Like the traditional flux map, processed Husimi flows can be integrated over lines and surfaces to reveal net current flow.

The processed Husimi formalism is able to identify several classical paths whose associated quantum waves are interfering: something that ordinary quantum flux cannot capture in any obvious way. This can be done if the dominant waves are sufficiently separated in k-space, i.e., if the momentum uncertainty of the coherent state can resolve between the trajectories, they can be found using the optimization (multi-modal analysis) algorithm.

The Husimi vectors in fig. 1 reveal that even a single plane wave can produce a range of Husimi vectors. Figure 2 demonstrates the multi-modal analysis (MMA) algorithm (see table 1) on the pure momentum state $\Psi _{{\rm C}}({\mathbf { r}}) = \alpha \cos \left ({\mathbf { k}}_{1}\cdot {\mathbf { r}}\right )+\beta \cos \left ({\mathbf { k}}_{2}\cdot {\mathbf { r}}\right )$, where the wave vectors are indicated by the white arrow overlay with k₁ pointing towards the upper-right and k₂ in the upper-left direction and α = 1 and β = 0.8. The Husimi vectors shown in fig. 2(a) have an uncertainty of Δk/k = 30%.

Zoom In Zoom Out Reset image size

Table 1:. Multi-modal analysis (MMA) algorithm.

1)	A set of templates on N wave vectors $\left \{ {\mathbf { k}}_{j}\right \}$ is created for the wave functions Ψ=eiktesti·r generated by the M wave vectors $\left \{ {\mathbf { k}}_{i}^{{\rm test}}\right \}$. Both sets of wave vectors lie along the dispersion contour. Each template is stored as a vector of values ui (length M) and corresponds to eq. (12) along kj.
2)	The metric given by di=v·ui, where the vector v represents the Husimi vector, is computed for each template in step 1).
3)	The maximum of the set $\left \{ d_{i}\right \}$ and the corresponding wave vector ktesti are stored.
4)	The contribution of the trajectory with wave vector ktesti is determined by the vector ${\mathbf { u}}_{i}\frac {d_{i}}{{\mathbf { u}}_{i}\cdot {\mathbf { u}}_{i}}$.
5)	The re-weighted template vector is subtracted from the Husimi vector, i.e., ${\mathbf { v}}\rightarrow {\mathbf { v}}-{\mathbf { u}}_{i}\frac {d_{i}}{{\mathbf { u}}_{i}\cdot {\mathbf { u}}_{i}}$.
6)	All negative elements of v are set to zero.
7)	Steps 1)–6) are repeated until the metric di dips below a threshold.
8)	The set of vectors $\left \{ d_{i}{\mathbf { k}}_{i}^{{\rm test}}\right \}$ are used to approximate the processed Husimi flow.

Although the calculations reported here are for very simple models, we believe that the technique is generic and can be applied to virtually any quantum wave function and model amenable to semiclassical analysis [14–16]. We begin by demonstrating the processed Husimi technique on a circular billiard, which, due to time-reversal symmetry, has zero flux. The Schrödinger equation for this system can be written in radial form as

$$\begin{eqnarray*} &&\frac{{\rm d}^{2}R(r)}{{\rm d}r^{2}}+\frac{1}{r}\frac{{\rm d}R(r)}{{\rm d}r}+\left(k^{2}-\frac{m^{2}}{r^{2}}\right)R(r)=0. \end{eqnarray*}$$

Solutions to this equation are simultaneous eigenstates of energy and angular momentum, and thus possess the good quantum numbers n (number of nodes in the radial direction) and m (number of angular nodes). Figure 3(b) shows one such state with n ≈ m, which corresponds to classical paths that bounce off the boundary at a consistent 30° angle (see schematic in fig. 3(d).

Zoom In Zoom Out Reset image size

Because the flux is uniformly zero, it is no help in understanding the dynamics of this system. In Fig. 3(c) we present the processed Husimi flow of this wave function by sampling coherent-state projections on a regularly spaced grid across the system. With the solid line, we show that these projections align perfectly with one of the classical paths indicated in fig. 3(d). In addition, each point in the Husimi map contains an additional set of Husimi vectors, which do not align with the path. Given that any state of a circularly symmetric system must correspond with infinitely many classical trajectories related by rotation, the processed Husimi at a particular point must reflect all rotated paths that intersect there. We indicate one classical path and its rotated version in grey in fig. 3(d). The "cross-hatching" pattern arises because two rotated classical paths intersect at any point, explaining the similar cross-hatching nodal patterns in the wave function.

Magnetic systems can also be analyzed using processed Husimi projections, thus shedding light on recent work examining flux vortices in quantum dots [17, 18]. Time-reversal symmetry in the circular billiard is broken by the magnetic field. To properly represent these states, both the momentum operator in the flux operator (eq. (1)) and the momentum term ik₀·r₀ in the coherent state (eq. (9)) must be modified to reflect the canonical transformation ${\mathbf { p}}\rightarrow {\mathbf { p}}-q{\mathbf { A}}/c$, where A is the magnetic potential.

For strong enough magnetic fields, classical trajectories are circular with radii corresponding to the cyclotron radius. However, the dynamics due to the cyclotron radius are not obvious from the wave function in fig. 4(b). Worse still, the flux map in fig. 4(a) (obtained by sampling at regularly spaced points in a grid similar to the one used for the processed Husimi flow) appears to indicate two circular paths but with strongly different radii.

Zoom In Zoom Out Reset image size

The processed Husimi flow in fig. 4(c) dramatically clarifies the dynamics by indicating classical paths consistent with a cyclotron radius approximately one-third of the radius of the system. The cyclotron orbits are offset from the center of the system; due to circular symmetry, there exist infinitely many such cyclotron orbits rotated around the system center. We show a subset of these orbits in the schematic in fig. 4(d), which closely parallels the processed Husimi flow and exhibits the correct cyclotron radius.

Why does the flux map fail to show the dominant classical paths? In fig. 5, we provide magnified views from the full set of coherent-state projections, the flux map, and the processed Husimi flow corresponding to the circles in fig. 4. We can model this point in the wave function according to the pure momentum state

$$\begin{equation} \Psi\left({\bf r}\right)=e^{i{\bf k}_{1}\cdot{\bf r}}+e^{i{\bf k}_{2}\cdot{\bf r}}, \end{equation} \tag{ 14 }$$

where k₁ and k₂ are indicated by the arrow overlays in fig. 5(c). Because of the uncertainty of each coherent-state projection, the full set of projections exhibits a finite spread around the generating wave vectors. The processed Husimi technique (red) retrieves two independent trajectories at this point, while the flux (blue) averages over them. The flux map amounts to summing the Husimi vectors, giving a total flow at each point that does not always correspond to the semiclassical dynamics of the system. With the processed Husimi flow, we now have a complementary representation for revealing the flow structure present in the wave function. Further, this formalism is not limited to quantum wave functions as there might be applications to optical waves. The application of processed Husimi techniques in optical currents [4] would be a significant step towards resolving a discordance that goes back to Young's interference experiments: the flux of a superposition of wave functions is not equal to the superposition of fluxes [19, 20].

Zoom In Zoom Out Reset image size

Processed Husimi maps also have implications for experiments measured in a fashion similar to angle-resolved photoemission spectroscopy (for a review, see [21]). In the ARPES setup, a focused photon beam on a sample emits electrons from the valence band. The energy of the photoemitted electrons incorporates both their bonding energies, which can be averaged over, and their kinetic energy, which depends on the angle of the beam with respect to the sample surface.

The ARPES response function behaves similarly to coherent-state projections with k₀ proportional to the beam angle. By rotating the beam angle around the same point of intersection, the response in different directions provides the momentum distribution of the wave function at that point. Perturbations from the known dispersion relation can then be inserted into eq. (13) to obtain the flux expectation value. While a narrow beam could measure the flux vector at the intersection point, it will be difficult to distinguish the occasional large perturbation measurements from noise. However, wider beams could capture additional terms from the Taylor expansion of the coherent state in eq. (10), producing more reliable measurements. Applying the technique at many points across the sample would then provide the processed Husimi map.

We have provided a new interpretation of the flux operator from the perspective of its eigenstates while connecting them to coherent-state projections at the limit of infinitesimal spatial spread. Away from this limit, we can use coherent-state projections to provide a map of the classical dynamics underlying a quantum wave function, permitting us to describe flow even in stationary states with zero flux. In systems with magnetic fields, we have shown that the flux maps actually correspond to the aggregate of such classical flows, which we are able to retrieve from the processed Husimi projections.

This research was conducted with funding from the Department of Energy Computer Science Graduate Fellowship program under Contract No. DE-FG02-97ER25308. MFB and EJH were supported by the Department of Energy, office of basic science (grant DE-FG02-08ER46513).