# Thermodynamics of Harmonic Bosons

Class Session 07 : Thermodynamics of Harmonic Bosons

# Introduction

Bosons are indistinguishable particles, which implies that the wave function describing an ensemble of bosons is symmetric with respect to any permutation of the individual positions of the bosons. Consequently, the density matrix of N non-interacting bosons does not simply factorise into a product of N single particle density matrices.

In this Class Session we examine some properties of the partition function of non-interacting bosonic particles.

# The partition function of N non-interacting bosons

## Canonical partition function of bosonic particles

• Write the density matrix of N distinguishable particles described by their positions (x1,...,xN) in terms of the single particle density matrix.
• Explain why the same relation does not hold for bosonic particles.
• Find a relation between the density matrix of N bosonic particles and the density matrix of N distinguishable particles.
• Let's now consider the canonical partition function ZN of N bosonic (i.e. indistinguishable) particles described by positions (x1,...,xN). Show that ZN writes as:

$Z_N= \frac 1{N!}\sum_P P\!Z_N^{\text{dist}}$

• where the sum is over the N! permutations P of (1,...,N), and PZNdist is the the partition function of N distinguishable particles where the final positions has been permuted by P:

$P\!Z_N^{\text{dist}} = \int d^N x\ \rho^{\text{dist}}(\{x_1,\ldots,x_N\},\{x_{P(1)},\ldots,x_{P(N)}\},\beta)$

• In this expression, ρdist is the density matrix of N distinguishable particles at inverse temperature is β = 1 / T.
• We note that the partition function writes ZN writes as ZN = YN / N! with YN a weight which writes as

$Y_N = \sum_P \text{weight}(P)$

• Any permutation P of (1,...,N) can be decomposed as a unique product of commuting cycle permutations. By partitionning the sum ΣP over the length k of the cycle to which the last element of the permutation P belongs, show the following recursion relation:

$Y_N = \sum_{k=1}^N z_k \binom{N-1}{k-1} (k-1)!\: Y_{N-k} \qquad \text{with}\quad Y_0=1$

• where zk is the single particle partition function at inverse temperature .
• Show now [Borrmann and Franke 1993] the following recursion relation for the partition ZN:

$Z_N= \frac 1N \sum_{k=1}^N z_k Z_{N-k} \qquad\qquad \qquad\qquad \qquad\qquad \qquad(1)$

• with Z0 = 1.
• In other words, the probability πk that the last particle belongs to a cycle of length k in the permutation P is:

$\pi_k = \frac 1N \frac {z_k}{Z_{N}}{Z_{N-k}} \qquad\qquad \qquad\qquad \qquad\qquad \qquad(2)$

## The case of bosons in an harmonic trap

We chose a harmonic trap the such that the energy levels in each of the three spatial directions of the trap are En = n. The (three-dimensional) single particle partition zk writes:

$z_k= \Big(\frac{1}{1-e^{-k\beta}}\Big)^3$

$T_\star = \frac{T}{N^{1/3}}$

## Algorithm

from pylab import *
from math import sqrt, sinh, tanh, exp
from random import uniform as ran, gauss
def z(beta,k):
sum = 1 / (1 - exp(-k*beta))**3
return sum
def canonic_recursion(beta,N):
Z = [1.]
for M in range(1, N+1):
Z.append( sum(Z[k] * z(beta,M-k) for k in range(M)) / M )
return Z
def pi_list_make(Z,M):
pi_list = [0] + [ z(beta,k)*Z[M-k]/Z[M]/M for k  in range(1,M+1) ]
pi_sum = [0]
for k in range(1, M+1):
pi_sum.append(pi_sum[k-1] + pi_list[k])
return pi_sum
def tower_sample(data,Upsilon): #naive tower sampling, cf. SMAC Sect. 1.2.3
for k in range(len(data)):
if Upsilon < data[k]: break
return k
def levy_harmonic_path(Del_tau,N): #
beta = N * Del_tau
xN = gauss(0.,1./sqrt(2*tanh(beta/2.)))
x = [xN]
for k in range(1,N):
Upsilon_1 = 1./tanh(Del_tau) + 1./tanh((N-k)*Del_tau)
Upsilon_2 = x[k-1]/sinh(Del_tau) + xN/sinh((N-k)*Del_tau)
x_mean = Upsilon_2 / Upsilon_1
sigma = 1. / sqrt(Upsilon_1)
x.append(gauss(x_mean,sigma))
return x

N = 512
T_star = .1
T = T_star * N**(1./3.)
beta = 1. / T
Z = canonic_recursion(beta,N)
M = N
x_config = []
y_config = []
while M > 0:
pi_sum = pi_list_make(Z,M)
Upsilon = ran(0,1)
k = tower_sample(pi_sum,Upsilon)
x_config += levy_harmonic_path(beta,k)
y_config += levy_harmonic_path(beta,k)
M -= k
figure(1)
axis('scaled')
xlim(-5.,5.)
ylim(-5.,5.)
plot(x_config,y_config,'ro')
xlabel('$x$')
ylabel('$y$')
title('3d bosons in trap (2d projection):'+' $N$= '+str(N)+' T* =' + str(T_star))
show()
• Familiarize yourself with the above algorithm, and run it for values of temperatures of T٭. What is the advantage of defining T٭?
• Modify the program such that, together with the scatter plot of particles, it also draws the longest cycle in the configuration. Plot snapshots at various values of T٭.
• Modify the program so that it computes the cycle-length distribution at the values of T٭ you studied before. Compare the values of a Monte Carlo sampling, and a direct plot of the distribution using equation (1) above. Show that the cycle distribution has a flat part, and comment.
• Textbooks of theoretical physics tell us that the quantum nature of bosonic fluids becomes noticeable when the inter-particle distance is of the order of the thermal de Broglie wavelength. We would rather say that the quantum nature of bosonic fluids becomes noticeable when long cycles appear. Why are the two statements closely related?

## Modified Algorithm

The following program enables to compare the theoretical and numerical distribution of cycles:
from pylab import *
from math import sqrt, sinh, tanh, exp
from random import uniform as ran, gauss
from matplotlib import pyplot
def z(beta,k):
sum = 1. / (1. - exp(-k*beta))**3
return sum
def canonic_recursion(beta,N):
Z = [1.]
for M in range(1,N+1):
Z.append(sum(Z[k] * z(beta,M-k) for k in range(M))/M)
return Z
def pi_list_make(Z,M,beta):
pi_list=[0] + [z(beta,k)*Z[M-k]*1.0/Z[M]/M for k  in range(1,M+1)]
pi_sum=[0]
for k in range(1,M+1):
pi_sum.append(pi_sum[k-1]+pi_list[k])
return pi_sum
def tower_sample(data,Upsilon):
for k in range(len(data)):
if Upsilon<data[k]: break
return k
#-------------------------------------------------------
# -- main program starts here, don't choose N too large
#-------------------------------------------------------
def run(T_star,N):
T=T_star*N**(1./3.)
beta=1./T
Z=canonic_recursion(beta,N)
M=N
lengths = []
while M > 0:
pi_sum=pi_list_make(Z,M,beta)
Upsilon=ran(0,1)
k=tower_sample(pi_sum,Upsilon)
lengths.append(k)
M-=k
i = lengths.index(max(lengths))
return lengths

N=64
Tstar=.3
M=4000
confs=[]
for m in range(M):
lengths=run(Tstar,N)
confs.append(lengths[0])
pyplot.hist(confs,N,range=(0,N+1),normed=True,histtype='step')
T=Tstar*N**(1/3.)
beta=1./T
Z=canonic_recursion(beta,N)
plot(range(1,N),[1./N*z(beta,k)*Z[N-k]/Z[N] for k in range(1,N)])
savefig("cycle_distribution_harmonic-bosons.png")
show()


 Figure 1. Cycle distribution at N=64: theortical (full line) and numerical (histogram) evaluations (M=4000 samples). Reduced temperature is 0.3

 Figure 2. Cycle distribution at N=64: theortical (full line) and numerical (histogram) evaluations (M=4000 samples). Reduced temperature is 0.5

 Figure 3. Cycle distribution at N=64: theortical (full line) and numerical (histogram) evaluations (M=4000 samples). Reduced temperature is 0.65