Class Session 06 : The Method of Images

Introduction

We are interested in the density matrix of a non-interacting particle living in a finite or semi-infinite interval of the real line, and subjected to various boundary conditions.

The free particle in a box

Consider a particle living in a box [0,L].

Case of the periodic box

  • We start from expression (3.11) of SMAC for the density matrix with periodic boundary conditions. Using the Poisson summation formula



  • show that



  • Check that



  • You can show this identity by direct computation and by using the evolution equation of the density matrix, together with the boundary conditions. Make the analogy with a particle diffusing on a ring.

Case of the box with hard walls (SMAC 3.1.3 pages 137-139)

  • What are the eigenvectors and eigenvalues of a free particle living inside a box [0,L]?
  • Compute the corresponding density matrix ρbox,L(x,x',β) at inverse temperature β.
  • Express ρbox,L(x,x',β) in terms of the free density matrix ρfree,L(x,x',β), using again the Poisson summation formula, as



  • This writes: ρbox, L(x,x',β) = ρper, 2L(x,x',β) − ρper, 2L(x, −x',β).

  • Explain why





Method of images for a single hard wall

Using a geometrical construction for the Feynman paths, it is possible to obtain the same result in few lines. We start by considering a simpler example: a particle living on the semi-infinite line [0,+[ with a unique hard wall at x = 0.
image_construction_wall.png
Figure 1. The free density matrix as a sum over two class of path.

  • Analysing figure 1 and thinking a little bit, show that one can write the density matrix of the particle at inverse temperature β as follows:



  • Compute the following integral and comment.



  • By using the same trick and decomposing paths in different classes (see figure 2 and 3 below), it is possible to rederive the expression of the density matrix in a box (see SMAC part 3.3.3 pages 155-157).

Method_of_images_and_the_particle_in_a_box.png
Figure 2. The free density matrix as a sum over three classes of paths.

density-matrix_left-right.png
Figure 3. Transformation from right to left class of paths.



Sampling in a box

We now sample the paths contributing to ρbox,L(x,x',β).
  1. First idea: start a Lévy construction and reject if the path is outside the box. Can you estimate the rejection rate?
  2. Smart Trick: it is better to sample the path on large scales first, and then to work one's way to small scale. Do you improve the acceptance rate?
  3. Imagine a rejection-free Monte-Carlo algorithm making use of this formula.

Idea number 1:
import random, pylab
import math
beta = 10
N = 500
Lbox = 10
dt = beta/float(N)
B=[0 for n in range(N+1)]
B[0] = 1.
B[N] = 3.
xlist = [n*dt for n in range(N+1)]
itot = 1
ireject = 0
while itot < N:
     DT = dt * (N-itot)
     mu = ( DT*B[itot-1] + dt*B[N] ) / ( dt+DT )
     sigma = 1. / math.sqrt( (dt+DT) / dt / DT)
     Bnew = random.gauss(mu, sigma)
     B[itot] = Bnew
     if Bnew < 0: itot,ireject = 0,ireject+1
     if Bnew > Lbox: itot,ireject = 0,ireject+1
     itot += 1
print beta, 'beta', ireject,'rejection number'
pylab.title('Sampling in a box')
pylab.plot(B,xlist, 'r-')
pylab.plot([0 for n in range(N+1)],xlist, 'k-')
pylab.plot([Lbox for n in range(N+1)],xlist, 'k-')
pylab.axis([-5, Lbox+5, 0, beta])
pylab.show()
exit()

Analogy with standard diffusion of a one dimensional Brownian particle

There is an interesting analogy between our quantum-mechanical non-interacting particle and the classical, stochastic, Brownian motion.
  1. Can you explicit this analogy? What is the signification of the spatial integrals of ρ(x,x',β) that you have computed?
  2. What are the boundary conditions conditions verified by the density matrix ρwall(x,x',β) given in equation (1)?
  3. Find the equivalent of equation (1) for reflecting boundary conditions.
from pylab import *
x0 = 4.
time = 500
steps = 200
xspace = linspace(0, 15, steps)
dt = 0.7
rhob = []
for n in range(time):
     rhob.append([exp(-(xspace[i]-x0)**2/2./dt)/sqrt(2*pi*dt)-exp(-(xspace[i]+x0)**2/2./dt)/sqrt(2*pi*dt) for i in range(len(xspace))])
     dt += 0.025
ion()#interactive mode
line, = plot(xspace, rhob[0])
plot(xspace, rhob[0])
plot([0,0],[1,0],'k-')
axis([xspace[0]-5, xspace[-1], 0, .5])
for n in range(len(rhob)):
    line.set_ydata(rhob[n])
    draw()


[Print this page]