Class Session 05 : The Evolution Operator

Introduction

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 β. In this Class Session, we examine the different problem of finding the time-evolution of a wave function Ψ(t) which starts from an initial state Ψ0 at time t=0, and evolves through the Schrödinger equation.



We study in this Class Session how to tackle this problem using the evolution operator.

A naive method: time discretization

Solving the Schrödinger evolution by remplacing the time derivative by a first order approximation



is simple to implement but raises problems (it is difficult to ensure the conservation of probability).

A better method: the time evolution operator


It is better to start from the operator solution of the Schrödinger equation



We observe that, putting formally β = 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 β.

Implementing the time evolution

We consider a particle of mass m=1, momentum p and position x, in a potential V(x). Its wave function Ψ(x,t) evolves with the Schrödinger equation of Hamiltonian



where p2 has the following representation in direct space




The Trotter formula

As in the equilibrium approach, the idea is to decompose the time interval t in N steps of duration Δt = t / N. The full operator of evolution can be written (without approximation) as a product of N operators



where each infinitesimal operator of evolution writes



In this expression, the potential part V(x) takes a simple form in real space, while the derivative part p2 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



To which order in Δt is that formula valid? Check this fact explicitely.

The algorithm

A possible scheme for the algorithm is thus:
  • Start from an initial wave function.
  • Discretise the time t in N steps of duration Δt = t / N .
  • At each time step, apply the Trotter formula by:
  1. in real space, multiplying by e −i Δt V(x) / 2,
  2. going to momentum space and multiply by eΔt p^2 / 2,
  3. going back to real space and multiply by e −i Δt V(x) / 2.


Notes for the implementation


  • The continuous space Fourier transforms write




In the algorithm, the space coordinate is also discretized and one uses the discrete Fourier transform.
  • An array may be used to encode the wave function: this simplifies the operations.
  • The discrete Fourier transform assumes the system to be periodic in space. This means, that the algorithms in the examples below actually simulate a wave function in a periodic potential landscape.

Examples

Oscillating in the harmonic potential

from pylab import *
steps = 400
space = linspace(-20,20,steps)
xspace = pspace = space       # real space and momentum space lattices
delta = space[1]-space[0]
T=40                          #final time
t_step=0.5
def fourierxp(phix):
    return array([integral(phix * exp(-1j * p * xspace)) for p in pspace])
def fourierpx(phip):
    return array([1./(2*pi)*integral(phip * exp(1j * x * pspace)) for x in xspace])
def integral(integrand):
    return sum(integrand)*delta
def time_step(psi0,pot):
    psi0 = exp(-1j * pot * t_step / 2) * psi0
    psi0 = fourierxp(psi0)
    psi0 = exp(-1j * pspace**2 * t_step / 2) * psi0
    psi0 = fourierpx(psi0)
    return exp(-1j * pot * t_step / 2) * psi0
def zeit_ent(psi0, pot):
    evaluation = psi0
    t = 0
    while t < T:
        print t
        t += t_step
        evaluation = concatenate((evaluation,time_step(evaluation[-1],pot)))
    return evaluation
pot=array([ [x**2/10 for x in xspace] ])
psi=array([ [2. * 1./sqrt(2*pi) * exp(-((x-1)/2)**2 / 2) for x in xspace] ])
psi = zeit_ent(psi, pot)
ion()
line, = plot(xspace, abs(psi[0])**2)
plot(xspace, pot[0]/4)
axis([space[0],space[-1],-.1,1])
for n in range(len(psi)):
    line.set_ydata(abs(psi[n])**2)
    draw()
show()

Harmonic-potential_wave-packet.png
An harmonic trap (in green) with the initial probability density (in blue).


We consider an harmonic (i.e. quadratic) potential with an inital Gaussian wave function.
  • What is your prediction for the time evolution of the probability density |Ψ(x,t)|2 ?
  • Check by coding the time evolution of the wave function in python.

The tunnel effect


tunnel-potential_wave-packet.png
The tunnel potential (in green) with the initial probability density (in blue).


  • What is your prediction for the time evolution of the probability density |Ψ(x,t)|2 ?
  • Check by coding the time evolution of the wave function in python.

Interferences in a Gaussian trap


gaussian-trap_wave-packet.png
The Gaussian trap potential (in green) with the initial probability density (in blue).


  • What is your prediction for the time evolution of the probability density |Ψ(x,t)|2 ?
  • Check by coding the time evolution of the wave function in python.


Example animations


Oscillating in the harmonic potential


Dev_ho.gif
The Gaussian trap potential (in green) with the time-dependent probability density (in blue).


The tunnel effect


Dev_tunnel.gif
The tunnel potential (in green) with the time-dependent probability density (in blue).


Interferences in a Gaussian trap


Maxwell_Potential.gif
The Gaussian trap potential (in green) with the time-dependent probability density (in blue).


Code version with gif-file for animation

from pylab import *
import time, os
steps = 400
space = linspace(-20,20,steps)
xspace = pspace = space       # real space and momentum space lattices
delta = space[1]-space[0]
T=40                         #final time
t_step=0.5
def fourierxp(phix):
    return array([integral(phix * exp(-1j * p * xspace)) for p in pspace])
def fourierpx(phip):
    return array([1./(2*pi)*integral(phip * exp(1j * x * pspace)) for x in xspace])
def integral(integrand):
    return sum(integrand)*delta
def time_step(psi0,pot):
    psi0 = exp(-1j * pot * t_step / 2) * psi0
    psi0 = fourierxp(psi0)
    psi0 = exp(-1j * pspace**2 * t_step / 2) * psi0
    psi0 = fourierpx(psi0)
    return exp(-1j * pot * t_step / 2) * psi0
def zeit_ent(psi0, pot):
    evaluation = psi0
    t = 0
    while t < T:
        #print t
        t += t_step
        evaluation = concatenate((evaluation,time_step(evaluation[-1],pot)))
    return evaluation
pot=array([ [x**2/10 for x in xspace] ])
psi=array([ [2. * 1./sqrt(2*pi) * exp(-((x-1)/2)**2 / 2) for x in xspace] ])
psi = zeit_ent(psi, pot)
line, = plot(xspace, abs(psi[0])**2)
plot(xspace, pot[0]/4)
axis([space[0],space[-1],-.1,1])
for n in range(len(psi)):
    line.set_ydata(abs(psi[n])**2)
    savefig("img%05d.png" % (n))
input = "img*.png"
output = "img.gif"
print("starting conversion")
os.system("convert -delay 50 -loop 0 %s %s" % (input, output))
os.system("rm %s" % input)
 


[Print this page]