We have seen in the lecture how to treat equilibrium quantum problems, using the density matrix of a system in contact with a thermal bath of inverse temperature $\beta$. In this Class Session, we examine the different problem of finding the time-evolution of a wave function $\psi(t)$ which starts from an initial state $\psi(0)$ at time $t=0$, and evolves through the Schrödinger equation.
$$i\frac{\partial \psi(t)}{\partial t} = H\psi(t)\quad(\text{with }\hbar=1)$$We study in this Class Session how to tackle this problem using the evolution operator.
Solving the Schrödinger evolution by remplacing the time derivative by a first order approximation
$$i\frac{\partial \psi(t)}{\partial t} \simeq \frac{\psi(t+\Delta t)−\psi(t)}{Δt}$$is simple to implement but raises problems (it is difficult to ensure the conservation of probability).
It is better to start from the operator solution of the Schrödinger equation
$$\psi(t)=\mathrm{exp}(−itH)\psi(0)$$We observe that, putting formally $\beta = i t$ (that is, working in imaginary time), one recovers the same evolution as that of the density matrix for an equilibrium system at inverse temperature $\beta$.
We consider a particle of mass $m=1$, momentum $p$ and position $x$, in a potential $V(x)$. Its wave function $\psi(x,t)$ evolves with the Schrödinger equation of Hamiltonian
$$H=\frac 12 p^2 + V(x)$$where $p^2$ has the following representation in direct space
$$p^2=−\frac{\partial^2}{\partial x^2}.$$As in the equilibrium approach, the idea is to decompose the time interval $t$ in $N$ steps of duration $\Delta t = t / N$. The full operator of evolution can be written (without approximation) as a product of $N$ operators
$$ \mathrm{exp}(−itH)=\prod_{k=1}^N\mathrm{exp}(−i \Delta t H) $$where each infinitesimal operator of evolution writes
$$ \mathrm{exp}(−i \Delta t H) = \mathrm{exp}[−i\Delta t (p^2/2+V(x))]$$.
In this expression, the potential part $V(x)$ takes a simple form in real space, while the derivative part $p^2$ takes a simple form in the momentum space (that is, after a Fourier transform of the space coordinate). To separate both contributions, we use the Trotter formula
$$ \mathrm{exp}[−i\Delta t (p^2/2+V(x))]\simeq \mathrm{exp}[−iV(x)\Delta t/2] \mathrm{exp}[−i \Delta t p^2/2] \mathrm{exp}[−iV(x) \Delta t /2]. $$Questions
A possible scheme for the algorithm is thus:
In the algorithm, the space coordinate is also discretized and one uses the discrete Fourier transform.
Defining some meshes for the $x$-space and $k$-space.
from pylab import *
a, b, N = -10,10,2**9
x_grid = linspace(a,b,N)
delta_x = x_grid[1]-x_grid[0]
k0 = -5 # -5 for harmonic -18 for tunnel effect
k_grid = linspace(k0,-k0,N)
delta_k = k_grid[1]-k_grid[0]
Fourier transform, matrix way (could be improved using FFT) or on the fly :
isqrt2pi = 1/sqrt(2*pi)
FT = zeros(shape=(N,N),dtype=complex)
for m in range(len(k_grid)):
for n in range(len(x_grid)):
FT[m,n] = delta_x*isqrt2pi*exp(-1j * k_grid[m] * x_grid[n])
iFT = zeros(shape=(N,N),dtype=complex)
for m in range(len(k_grid)):
for n in range(len(x_grid)):
iFT[n,m] = delta_k*isqrt2pi*exp(+1j * k_grid[m] * x_grid[n])
def ft(phix):
#return dot(FT,phix)
return isqrt2pi * array([ delta_x * sum(phix * exp(-1j * k * x_grid)) for k in k_grid])
def ift(phik):
#return dot(iFT,phik)
return isqrt2pi * array([ delta_k * sum(phik * exp( 1j * x * k_grid)) for x in x_grid])
def norm_x(psi): return sum(delta_x*abs(psi)**2)
def norm_k(psi): return sum(delta_k*abs(psi)**2)
Some wave-functions and basic checks on Fourier transform
def gaussian_wave(x_,x0_=0,sigma=1,k0_=0):
return sqrt(isqrt2pi/sigma) * exp(-( (x_-x0_)/(2*sigma) )**2 + 1j*k0_*x_)
def plotwf(psi):
fig = figure()
ax1 = fig.add_subplot(211)
ax1.plot(x_grid, psi.real, 'r+-',label=r'$\Re\psi(x)$')
ax1.plot(x_grid, psi.imag, 'g+-',label=r'$\Im\psi(x)$')
ax1.plot(x_grid, abs(psi)**2, 'b-',label=r'$|\psi(x)|^2$')
line_x, = ax1.plot(x_grid, abs(psi)**2, 'b')
line_x.set_linewidth(2)
y_min = min(min(psi.real),min(psi.imag))
y_max = max(max(psi.real),max(psi.imag),max(abs(psi)**2))
ax1.axis([x_grid[0],x_grid[-1],y_min,y_max])
ax1.set_xlabel('$x$')
ax1.legend(prop=dict(size=16))
ax2 = fig.add_subplot(212)
line_k, = ax2.plot(k_grid, abs(ft(psi))**2, 'b', label=r'$|\psi(k)|^2$')
line_k.set_linewidth(3)
ax2.set_xlabel('$k$')
ax2.legend(prop=dict(size=16))
ax2.axis([k_grid[0],k_grid[-1],-.1,1])
gwp = gaussian_wave(x_grid,x0_=0,sigma=1,k0_=2)
print(norm_x(gwp))
# checking Parseval-Plancherel
igwp= ft(gwp)
print(norm_k(igwp))
# checking fourier back
print( norm_x(ift(igwp)-gwp) )
#plot(x_grid,abs(gwp)**2,'+-',x_grid,gwp.real,'r-')
#plot(k_grid,abs(igwp)**2,'+-')
plotwf(gwp)
show()
Time evolution simulation
def time_step(psi,pot,dt):
psi = exp(-1j * pot * dt / 2) * psi
psi = ft(psi)
psi = exp(-1j * k_grid**2 * dt / 2) * psi
psi = ift(psi)
return exp(-1j * pot * dt / 2) * psi
def time_evolution(psi0, pot, times):
results = [ psi0 ]
tp = times[0]
for t in times[1:]:
dt, tp = t-tp, t
results.append(time_step(results[-1],pot,dt))
return results
import matplotlib.animation as animation
def animate_psi_of_t(psi0, pot, times):
psi_of_t = time_evolution(psi0, pot, times)
fig = figure()
ax1 = fig.add_subplot(211)
ax1.plot(x_grid, pot, 'k')
ax1.plot(x_grid, abs(psi0)**2, 'b-')
line_x, = ax1.plot(x_grid, abs(psi0)**2, 'b',label=r'$|\psi(x)|^2$')
line_x.set_linewidth(3)
ax1.axis([x_grid[0],x_grid[-1],-.1,0.7])
ax1.set_xlabel('$x$')
ax1.legend(prop=dict(size=16))
ax2 = fig.add_subplot(212)
ax2.plot(k_grid, abs(ft(psi0))**2, 'b-')
line_k, = ax2.plot(k_grid, abs(ft(psi0))**2, 'b',label=r'$|\psi(k)|^2$')
line_k.set_linewidth(3)
ax2.axis([k_grid[0],k_grid[-1],-.1,1.2])
ax2.set_xlabel('$k$')
ax2.legend(prop=dict(size=16))
def animate(i):
line_x.set_ydata(abs(psi_of_t[i])**2)
line_k.set_ydata(abs(ft(psi_of_t[i]))**2)
return (line_x, line_k)
ani = animation.FuncAnimation(fig, animate, np.arange(len(psi_of_t)), interval=50)
show()
omega = sqrt(2/10)
def potential(x):
return 0.5*(omega*x)**2
def coherent_state(x_,alpha=1.0):
x0_ = sqrt(2/omega)*alpha.real
sigma = 1/sqrt(2*omega)
k0_ = sqrt(2*omega)*alpha.imag
return sqrt(isqrt2pi/sigma) * exp(-( (x_-x0_)/(2*sigma) )**2 + 1j*k0_*x_)
pot = potential(x_grid)
psi0 = gaussian_wave(x_grid,x0_=0,sigma=.8,k0_=.6)
#psi0 = coherent_state(x_grid,alpha=1.0)
times = arange(0,20,0.15)
animate_psi_of_t(psi0, pot, times)
We consider an harmonic (i.e. quadratic) potential with an inital Gaussian wave function.
def barrier(x,a=0,b=0.5,height=40):
if a<=x<=b: return height
else: return 0
pot = array([ barrier(x) for x in x_grid ])
psi0 = gaussian_wave(x_grid,x0_=-3,sigma=.4,k0_=10)
times = arange(0,2,0.005)
animate_psi_of_t(psi0, pot, times)