The Method of Images
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

  • 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. You will obtain:

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

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

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

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?
Figure 4. Path construction in a box, on increasingly finer scales

The algorithm
import random, pylab
import math
%matplotlib inline
beta = 10
N = 500
Lbox = 10
dt = beta/float(N)
B=[0 for n in range(N+1)]
B[0] = 1.
B[N] = 5.
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])

  1. Imagine a rejection-free Monte-Carlo algorithm making use of this formula, for example for position x4 in Fig. 4 above:

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 for the Brownian classical particles.

from pylab import *
x0 = 4.
time = 2000
steps = 800
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])
axis([xspace[0]-5, xspace[-1], 0, .5])
for n in range(len(rhob)):

[Print this page]