Solving Schrodinger’s equation: 
Coherent states

In this exercise, we’ll be exploring coherent states -- a solution to the quantum harmonic oscillator of the evolution of a displaced ground-state wavefunction. Again, we start by defining the constants and the ground state wavefunction. We’ll use McEuen’s buckyballs: m = 60 (Subscript[m, C]) = 60 * 12 (Subscript[m, p]) in a harmonic well with frequency one THz:

ℏ = ;
ω = ;
protonMass = ;
m = ;
a0= ;(* RMS width of Gaussian *)

We define ψ on a lattice of Np points x, from -L/2 to L/2, and Nt times separated by dt = Period / 100.

L = 30 a0;
Np = ;
dx = ;
x = ;
Period = ;
dt = ;
T = 2 Period;
Nt = Round[...]; (* Avoid rounding errors *)
times = Table[m dt, {m, 0, Nt}];

We solve Schrodinger’s equation by using the Baker–Campbell–Hausdorff formula to approximate the time evolution operator, applying the potential energy for dt/2, then the kinetic energy by dt, and then the potential energy by dt/2:
U(dt) = Upot(dt/2) Ukin(dt) Upot(dt/2)
 (e^(- i (p^2/2m + V(x)) dt/ℏ)) ~ (e^(- i V(x)dt /2ℏ)) (e^(-i (p^2/2m) dt/ℏ)) (e^(- i V(x)dt /2ℏ))
If ψ[x] is expressed in real space, applying the potential energy is just multiplication in real space:

V = ; 
UpotDtOver2 = ;

If ψ[k] is expressed in Fourier space, applying the kinetic energy is multiplication, as we saw for the free particle. Again, the FFT returns ψ[k] at points separated by dk = 2π/L: 
     k= {0, dk, 2 dk, ..., (Np/2) k, -(Np/2 -1) k, ..., -2 dk, -dk}

dk = ;
k = Join[];
k2 = ; 
UkinTildeDt = ;

Now we define the initial wavefunction ψ[0]. (The zero tells us this is the zero’th time-step.) . Let’s start in the ground state:

ψ[0] = ;
ψMax = Max[Abs[ψ[0]]];

We make a useful plotting function, that plots ψ[n] (real part, imaginary part, and absolute value). (Why don’t we plot the probability density |(ψ^2)(x,t)| along with the real and imaginary parts?)

Plotψ[n_] := ListPlot[{Transpose[{x, Re[ψ[n]]}], Transpose[{x, ...}], Transpose[{x, ...}]}, Joined->True, PlotRange->{- 1.1 ψMax,1.1ψMax}]

Plotψ[0]

We now evolve in small time steps dt, storing 
ψ[n](x) = ψ(x, n dt).

For[n = 1, n < Nt, n++,
ψ[n] = ... InverseFourier[...* Fourier[... ψ[n - 1]]];]

Plot the real part, the imaginary part, and the absolute value of ψ at time t = Period/5, hence n = Nt/10. (Why don' t we plot the probability density |ψ(x)((|)^2) along with the real and imaginary parts?)

Plotψ[...]

We use ListAnimate to see the evolution.

ListAnimate[Table[Plotψ[n],{n,0,Nt-1}]]

Note that the probability distribution doesn’t evolve with time; the ground state is a stationary state. What happens to the real and imaginary parts? Why?

Now let’s evolve from an initial condition that is the ground state displaced by 10 a0:

x0 = ...;
ψ[0] = ...; 
For[...];

Again, plot the evolution after time t=Period/5

Plotψ[...]

How does the probability density evolve, compared to a classical mass in the same harmonic oscillator?

ListAnimate[...]

Translating the ground state of the harmonic oscillator, either in position or in momentum, yields a coherent state. Coherent states are great tools for harmonic systems (like optics). They travel in phase space like classical particles. One can show that they are eigenstates of the annihilation operator. My graduate quantum course was filled with them (perhaps because the instructor, Roy Glauber, got the Nobel prize for developing them).