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:



  • 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:



  • 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



  • 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:



  • 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:



  • 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:





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:



We define a reduced temperature T٭ as (see also SMAC p.202):



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

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

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

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


[Print this page]

References

  • Borrmann P., Franke G. (1993) Recursion formulas for quantum statistical partition functions, Journal of Chemical Physics 98, 2484–2485
  • Landsberg P. T. (1961) Thermodynamics with quantum statistical illustrations, Interscience Publishers
  • SMAC part 4.2.3 - 4.2.6