Tutorial 05 : The Evolution Operator

Solving the time-dependent Schrödinger equation

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.

A naive method: time discretization

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

A better method: the time evolution operator

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$.

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 $\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}.$$

The Trotter formula

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

  1. To which order in $\Delta 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 $\Delta t = t / N$.
  • At each time step, apply the Trotter formula by:
    1. in real space, multiplying by $e^{ −i \Delta t V(x) / 2}$,
    2. going to momentum space and multiply by $e^{ −i \Delta t p^2 / 2}$,
    3. going back to real space and multiply by $e^{ −i \Delta t V(x) / 2}$.

Notes for the implementation

  • The continuous space Fourier transforms write $$ \psi(p)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}dx\,\mathrm{e}^{−ipx}\psi(x) \quad\text{ and }\quad \psi(x)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}dp\,\mathrm{e}^{ipx} \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

Defining some meshes for the $x$-space and $k$-space.

In [ ]:
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 :

In [ ]:
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])
In [ ]:
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

In [ ]:
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

In [ ]:
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()
In [ ]:
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.

  1. code the time evolution of the wave function in python.
  2. try to start from a coherent state. $$\psi_{\alpha}(x)=\left(\frac{m\omega}{\pi\hbar}\right)^{1/4}e^{-\frac{m\omega}{2\hbar}\left(x-\sqrt{\frac{2\hbar}{m\omega}}\Re[\alpha]\right)^2+i\sqrt{\frac{2m\omega}{\hbar}}\Im[\alpha]x} $$

Tunnel effet

In [ ]:
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)
In [ ]: