Animation
Physical Interpretation
The initial wavefunction is relatively widely spread in position space and has large support in the regions away from the origin, where the potential is larger. As time evolves, the distribution squeezes towards the origin, where the potential is minimized. This is the quantum analog of a classical particle that is attracted to the minimum of the potential due to a force F=-dV/dx. As the spatial distribution squeezes towards the centre x=0, its variance \sigma_x^2 becomes smaller. That is, the location of the particle becomes better defined. On account of Heisenberg’s uncertainty principle, however, this relativily well-defined position implies a large spread \sigma_p^2 in the momentum distribution. Subsequently, the large spread in momentum implies that the uncertainty in the position will be larger after the next time interval and the sequence repeats again: this results in the oscillatory pattern over time.
Note that, in contrast to the animation on free Decaying Gaussian Wavepacket, in this case the momentum spread is time-dependent as well. This is due to the fact that we have a spatially varying potential in this case, classically corresponding to a force.
Theoretical Context
We consider an initial state of the form:
\begin{equation*}
\Psi(x,0)=A\;x^2 e^{-x^2/2},
\end{equation*}trapped inside an harmonic potential of the form V(x)=x^2/2. For convenience, we work in natural units where \hbar=m=\omega\equiv 1. We wish to find the time-evolution of this initial state. In principle, this amounts to solving the Schrödinger equation, with the initial state as boundary condition, which is not a trivial task.
Solution in Real Space
Fortunately, the initial state above can be written as a superposition of energy eigenstates of the harmonic oscillator \psi_n(x). Once the correct superposition is found, the time-evolution is trivially obtained by appending the appropriate factors of e^{-i\omega_n t}, where \omega_n=n+1/2 is the energy of the n-th eigenstate. The eigenstates of the harmonic oscillator are the Hermite functions:
\begin{equation*}
\psi_n(x)=A_n\;H_n(x)\;e^{-x^2/2},
\end{equation*}where H_n(x) are the Hermite polynomials and A_n\equiv 1/(\sqrt{\pi}n!2^n)^{1/2} is the normalization factor. The first few Hermite polynomials are:
H_0(x)=1,\quad\quad\quad H_1(x)=2x,\quad\quad\quad H_2(x)=4x^2-2.
A general superposition \Psi(x,t) has the form:
\begin{equation*}
\Psi(x,t)=\sum_n c_n\;\psi_n(x)\;e^{-i\omega_n t}.
\end{equation*}Note that our initial wavefunction \Psi(x,0) can be written as a linear combination of the 0-th eigenstate (i.e. ground state) and the 2-nd eigenstate (i.e. second excited state), provided that we choose the weights c_{0,2} appropriately. To get a wave function proportional to x^2e^{-x^2/2}, the contribution of the ground state must cancel the part of the 2-nd eigenstate that is not proportional to x^2. That is, we must solve:
\begin{equation*}
c_0A_0 - 2c_2A_2 = 0,
\end{equation*}which gives c_2=\sqrt{2}c_0. Combining this with the constraint that c_0^2+c_2^2=1 gives us c_0=\sqrt{1/3} and c_2=\sqrt{2/3}. The wavefunction at any time t>0 is given by:
\begin{equation*}
\Psi(x,t)=c_0\;\psi_0(x)\;e^{-i\omega_0t}+c_2\;\psi_2(x)\;e^{-i\omega_2t},
\end{equation*}where you can check explicitly that we recover the initial form of the wavefunction at t=0 for the values found for the coefficients. Since we have written the wavefunction in terms of energy eigenstates, the expectation value for the energy is \langle E\rangle=c_0^2\omega_0+c_2^2\omega_2=11/6.
Solution in Momentum Space
So far, we found the solution in real space. To find the solution in momentum space, we have to Fourier transform the solution \Psi(x,t). Note that in the Fourier transform, position x and wavenumber k are conjugate variables, but since we have set \hbar=1, we know that momentum is identical to wavenumber (recall p=\hbar k). Hence momentum and position become conjugate variables. We denote the Fourier transform of f(x) as \tilde{f}(p)\equiv \mathcal{F}[f(x)](p) and use the convention:
\tilde{f}(p)=\int\frac{dp}{\sqrt{2\pi}}\;e^{-ipx}\;f(x). To find the transform of \Psi(x,t), note that our task is essentially to find the transform of the Hermite functions \psi_n(x) for the cases n=0,2, since all spatial dependence in contained in them. The Fourier transform of the Hermite functions is simply given by:
\begin{equation*}
\tilde{\psi}_n(p)=(-i)^n\psi_n(p).
\end{equation*}The Fourier transform of \Psi(x,t) then simply becomes:
\begin{equation*}
\tilde{\Psi}(p,t)=c_0\;\psi_0(p)\;e^{-i\omega_0t}-c_2\;\psi_2(p)\;e^{-i\omega_2t}.
\end{equation*}Properties of the Distributions
The probability density is given by |\Psi(x,t)|^2 and |\tilde{\Psi(x,t)}|^2 in real and momentum space, respectively. In terms of the Hermite functions, we have:
\begin{align*}
|\Psi(x,t)|^2 &= c_0^2\;|\psi_0(x)|^2+c_2^2\;|\psi_2(x)|^2+2c_0c_2\;\psi_0(x)\psi_2(x)\cos((\omega_2-\omega_0)t), \\
|\tilde{\Psi}(p,t)|^2 &= c_0^2\;|\psi_0(p)|^2+c_2^2\;|\psi_2(p)|^2-2c_0c_2\;\psi_0(p)\psi_2(p)\cos((\omega_2-\omega_0)t).
\end{align*}By inspection of the animations, the distributions are symmetric at all times, so that \langle x\rangle=\langle p\rangle=0. The variance is then simply \sigma_p^2=\langle p^2\rangle and \sigma_x^2=\langle x^2\rangle, which can be written as:
\begin{align*}
\sigma_x^2(t) & = c_0^2\;I_{0,0}^{(x)}+c_2^2\;I_{2,2}^{(x)}+2c_0c_2\;I_{0,2}^{(x)}\cos((\omega_2-\omega_0)t), \\
\sigma_p^2(t) & = c_0^2\;I_{0,0}^{(p)}+c_2^2\;I_{2,2}^{(p)}-2c_0c_2\;I_{0,2}^{(p)}\cos((\omega_2-\omega_0)t),
\end{align*}where the integral I_{m,n}^{(z)}, where z is the integration variable, is given by:
\begin{equation*}
I_{m,n}^{(z)}\equiv \int_{-\infty}^{+\infty}dz\;\psi_m(z) \;z^2\;\psi_n(z).
\end{equation*}For a closed form solution to this integral, see e.g. this link. We only require I_{0,0}^{(z)}=1/2,\;I_{2,2}^{(z)}=5/2 and I_{2,2}^{(z)}=5/2. Substituting the values for the coefficients as well, we find:
\begin{equation*}
\sigma_x^2(t)=\frac{11}{6}+\frac{2}{3}\cos((\omega_2-\omega_0)t),\quad\quad\quad \sigma_p^2(t)=\frac{11}{6}-\frac{2}{3}\cos((\omega_2-\omega_0)t).
\end{equation*}