In this lecture, we treat more advanced aspects of quantum mechanics: Bosons and Bose-Einstein condensation using the density matrix formalism. Once more the path integral representation of the density matrix will allow us to design new and efficient Monte Carlo algorithms.

Program:

From a single quantum particle to many quantum particles : the bosons and the Bose-Einstein condensation

Distinguishable and indistinguishable non-interacting particles: Energy levels and wave functions

Density matrix for many distinguishable non-interacting particles

Density matrix for bosons

Path integrals for bosons


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