In this lecture, we treat more advanced aspects of quantum mechanics, using the path integral representation: Interacting quantum particles, ideal bosons, and even interacting bosons. In all our examples, we stress the relationship of the Monte Carlo framework between quantum and classical problems.

Variations on the Free Path integral


Path correlations


The trivial Lévy program

see SMAC section 3.5.1

Fourier representation


Partition functions and density matrices as "sums of paths"


see SMAC section 3.4.1

Interacting particles


see SMAC section 3.4.1

Path integrals for bosons


To derive the path integral for bosons, we again use the position representation of operators and wave functions. We must simply remark that the density matrix is given by

Given a quantum system with an orthonormal set of wave functions

and with energies Ek, its density matrix at inverse temperature

is defined as


This remains true for bosons, if we use symmetric wave functions or for fermions, with anti-symmetric wave functions.

Markov-chain algorithm for bosons in a three-dimensional trap

import random, pylab
from math import sqrt, sinh, tanh, exp
 
def levy_harmonic_path(Del_tau, N):
    beta = N * Del_tau
    xN = tuple([random.gauss(0.0, 1.0 / sqrt(2.0 *
                  tanh(beta / 2.0))) for k in range(3)])
    x = [xN]
    for k in range(1,N):
        Upsilon_1 = 1.0 / tanh(Del_tau) + 1.0 / tanh((N - k) * Del_tau)
        Upsilon_2 = [x[k - 1][l] / sinh(Del_tau) + xN[l] /
                     sinh((N - k) * Del_tau) for l in range(3)]
        x_mean = [Upsilon_2[l] / Upsilon_1 for l in range(3)]
        sigma = 1.0 / sqrt(Upsilon_1)
        dummy = [random.gauss(x_mean[l], sigma) for l in range(3)]
        x.append(tuple(dummy))
    return x
 
def rho(x, xp, beta):
    Upsilon_1 = sum((x[l] + xp[l]) ** 2 / 4.0 *tanh(beta / 2.0) for l in range(3))
    Upsilon_2 = sum((x[l]-xp[l]) **2 / 4.0 / tanh(beta / 2.0) for l in range(3))
    return exp(- Upsilon_1 - Upsilon_2)
 
N = 512
T_star = 0.2
T = T_star * N ** (1.0 / 3.0)
beta = 1.0 / T
dd = {}
for k in range(N):
    a = levy_harmonic_path(beta, 1)
    dd[a[0]] = a[0]
for iter in range(30000):
    if iter%100 == 0: print iter
    perm=[]
    a = random.choice(dd.keys())
    while True:
        perm.append(a)
        b=dd.pop(a)
        if b == perm[0]: break
        else: a = b
    perm=levy_harmonic_path(beta,len(perm))
    dd[perm[-1]]=perm[0]
    for k in range(len(perm) - 1):
        dd[perm[k]] = perm[k + 1]
    a_1 = random.choice(dd.keys())
    b_1 = dd.pop(a_1)
    a_2 = random.choice(dd.keys())
    b_2 = dd.pop(a_2)
    weight_new = rho(a_1, b_2, beta) * rho(a_2, b_1, beta)
    weight_old = rho(a_1, b_1, beta) * rho(a_2, b_2, beta)
    if (random.random()) < weight_new / weight_old:
        dd[a_1] = b_2
        dd[a_2] = b_1
    else:
        dd[a_1] = b_1
        dd[a_2] = b_2
dummy = dd.keys()
x_config = [a[0] for a in dd.keys()]
y_config = [a[1] for a in dd.keys()]
pylab.figure(1)
pylab.axis('scaled')
pylab.xlim(-5.0, 5.0)
pylab.ylim(-5.0, 5.0)
pylab.plot(x_config, y_config, 'ro')
pylab.xlabel('$x$')
pylab.ylabel('$y$')
pylab.title('3d bosons in trap (MCQMC projection):'+' $N$= '+str(N)+' T* =' + str(T_star))
pylab.show()