# the evolution operator

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.

$i \frac{\partial}{\partial t} \Psi(t) = H \Psi(t) \quad\quad\text{(with \hbar=1)}$

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

$\frac{\partial}{\partial t} \Psi(t) \simeq \frac{\Psi(t+\Delta t) - \Psi(t)}{\Delta t}$

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

$\Psi(t) = \exp( - i t H) \Psi(0)$

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

$H=\frac{p^2}{2} + V(x)$

where p2 has the following representation in direct space

$p^2 = -\frac{\partial^2}{\partial x^2}\:.$

## 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

$\exp(- i t H) =\prod_{k=1}^N \exp( - i ~\Delta t~ H)$

where each infinitesimal operator of evolution writes

$\exp( - i ~\Delta t~ H) =\exp\Big[ - i ~\Delta t~ ( \frac{p^2}{2} + V(x)) \Big]\:.$

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

$\exp\Big[ - i ~\Delta t~ ( \frac{p^2}{2} + V(x)) \Big] \simeq \exp\Big[- i V(x)~\Delta t/2 \Big] ~ \exp\Big[-i ~\Delta t\frac{ p^2}{2}\Big] ~\exp\Big[ - i V(x)~\Delta t~/2\Big]\:.$

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

$\Psi(p) = \int_{-\infty}^{\infty} dx \exp(- i p x) \Psi(x)$
$\Psi(x) = \frac{1}{2 \pi} \int_{-\infty}^{\infty} dp \exp( i p x) \Psi(p)$

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()

 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

 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

 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

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

## The tunnel effect

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

## Interferences in a Gaussian trap

 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)