# 07_Quantum_Bosons

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:

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

```