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:

in real space, multiplying by e−i Δt V(x) / 2,

going to momentum space and multiply by e− Δtp^2 / 2,

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.5def fourierxp(phix):
returnarray([integral(phix * exp(-1j * p * xspace))for p in pspace])def fourierpx(phip):
returnarray([1./(2*pi)*integral(phip * exp(1j * x * pspace))for x in xspace])def integral(integrand):
returnsum(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 =0while t < T:
print t
t += t_step
evaluation = concatenate((evaluation,time_step(evaluation[-1],pot)))return evaluation
pot=array([[x**2/10for 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 inrange(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 *
importtime,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.5def fourierxp(phix):
returnarray([integral(phix * exp(-1j * p * xspace))for p in pspace])def fourierpx(phip):
returnarray([1./(2*pi)*integral(phip * exp(1j * x * pspace))for x in xspace])def integral(integrand):
returnsum(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 =0while t < T:
#print t
t += t_step
evaluation = concatenate((evaluation,time_step(evaluation[-1],pot)))return evaluation
pot=array([[x**2/10for 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 inrange(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)

Class Session 05 : The Evolution Operator## Table of Contents

## 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 timet=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 approximationis 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 inimaginary 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 massm=1, momentumpand positionx, in a potentialV(x). Its wave functionΨ(x,t) evolves with the Schrödinger equation of Hamiltonianwhere

p2 has the following representation in direct space## The Trotter formula

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

In this expression, the potential part

V(x) takes a simple form in real space, while the derivative partp2 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 formulaTo which order in

Δtis that formula valid? Check this fact explicitely.## The algorithm

A possible scheme for the algorithm is thus:tinNsteps of durationΔt=t/N .e−iΔt V(x) / 2,e−Δtp^2 / 2,e−iΔt V(x) / 2.## Notes for the implementation

In the algorithm, the space coordinate is also discretized and one uses the discrete Fourier transform.

## Examples

## Oscillating in the harmonic potential

We consider an harmonic (i.e. quadratic) potential with an inital Gaussian wave function.

Ψ(x,t)|2 ?## The tunnel effect

Ψ(x,t)|2 ?## Interferences in a Gaussian trap

Ψ(x,t)|2 ?## Example animations

## Oscillating in the harmonic potential

## The tunnel effect

## Interferences in a Gaussian trap

## Code version with gif-file for animation

[Print this page]