07_Quantum_Bosons

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=

The trivial Lévy program
see SMAC section 3.5.1

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 math \psi_k(x) math and with energies //E// k, its density matrix at inverse temperature math \beta=1/k_B T math is defined as math \rho(x,x',\beta)=\sum_k e^{-\beta E_k}\psi_k(x)\psi^*_k(x'). math

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
code format="python" 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

code