Calculations to accompany Section III of:¶
A computational introduction to quantum statistics using harmonically trapped particles¶
Martin Ligare, Department of Physics & Astronomy, Bucknell University¶
Abstract
In a 1997 paper Moore and Schroeder argued that the development of student understanding of thermal physics could be enhanced by computational exercises that highlight the link between the statistical definition of entropy and the second law of thermodynamics [Am. J. Phys. 65, 26 (1997)]. I introduce examples of similar computational exercises for systems in which the quantum statistics of identical particles plays an important role. I treat isolated systems of small numbers of particles confined in a common harmonic potential, and use a computer to enumerate all possible occupation-number configurations and multiplicities. The examples illustrate the effect of quantum statistics on the sharing of energy between weakly interacting subsystems, as well as the distribution of energy within subsystems. The examples also highlight the onset of Bose-Einstein condensation in small systems.
In [4]:
import numpy as np
from sympy.utilities import iterables
from scipy import special
import matplotlib as mpl
import matplotlib.pyplot as plt
In [5]:
# Following is an Ipython magic command that puts figures in notebook.
%matplotlib notebook
# M.L. modifications of matplotlib defaults
# Changes can also be put in matplotlibrc file,
# or effected using mpl.rcParams[]
mpl.style.use('classic')
plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('figure', autolayout = True) # Adjusts supblot parameters for new size
In [6]:
def mult_B(npa):
'''Calculate multiplicity for bosons in 3D trap with occupation-
number configuration in numpy array npa. See Eq. (13).'''
m = npa[0]
dm = (m+1)*(m+2)/2
nm = npa[1] # NOTE: broadcasting over configuration
a = special.binom(nm + dm -1, nm)
return np.prod(a)
def mult_D(npa):
'''Calculate multiplicity for distinguishable particle in 3D trap
with occupation-number configuration given in numpy array npa.
See Eq. (14).'''
m = npa[0]
dm = (m+1)*(m+2)/2
nm = npa[1]
n0 = n - np.sum(nm)
a = dm**nm/special.factorial(nm)
#print(n,dm,nm,n0,a)
return special.factorial(n)*np.prod(a)/special.factorial(n0)
Bosons¶
In [8]:
n = 20
q_max = 80
pop0_B = np.zeros(q_max+1)
pop0_B[0] = n
s_B = np.zeros(q_max+1)
for q in range(1,q_max+1):
ex0_B = 0
omega_B = 0
for p in iterables.partitions(q,n):
npa = np.array(list(p.items())).T
#print(npa)
omega1 = mult_B(npa)
ex0_B += omega1*(n - np.sum(npa[1]))
omega_B += omega1
pop0_B[q] = ex0_B/omega_B
s_B[q] = np.log(omega_B)
#print(q) # progress monitor
#print(ex0/omega, np.log(omega))
Distinguishable particles¶
In [9]:
n = 20
q_max = 80
pop0_D = np.zeros(q_max+1)
pop0_D[0] = n
s_D = np.zeros(q_max+1)
for q in range(1,q_max+1):
ex0_D = 0
omega_D = 0
for p in iterables.partitions(q,n):
npa = np.array(list(p.items())).T
#print(npa)
omega1 = mult_D(npa)
ex0_D += omega1*(n - np.sum(npa[1]))
omega_D += omega1
#print(p,omega1)
pop0_D[q] = ex0_D/omega_D
s_D[q] = np.log(omega_D)
#print(q) # progress monitor
Entropy¶
In [16]:
plt.figure()
x = np.linspace(0,q_max, q_max+1)
y1 = s_B
y2 = s_D
plt.xlim(0,80)
plt.ylim(0,100)
plt.scatter(x,y1, label='bosons')
plt.scatter(x, y2, label='distinguishable', color='red')
plt.xlabel('Energy (units: $\epsilon$)')
plt.ylabel('Entropy/$k$')
plt.title('Manuscipt Fig. 9')
plt.legend(loc='upper left');
Calculate temperatures¶
\begin{eqnarray*} \frac{1}{T_i} &=& \left.\frac{dS}{dq}\right\vert_i \\ &\simeq& \frac{S_{i+1} - S_{i-1}}{2\epsilon} \end{eqnarray*}Solving for $T_i$ gives $$ T_i = \frac{2\epsilon}{S_{i+1} - S_{i-1}}. $$
In [10]:
t_B = 2/(s_B[2:] - s_B[:len(s_B)-2])
t_D = 2/(s_D[2:] - s_D[:len(s_D)-2])
In [17]:
plt.figure()
x = np.linspace(0,q_max, q_max+1)
y1 = pop0_B
y2 = pop0_D
plt.scatter(x, y1)
plt.scatter(x, y2, color='red')
plt.ylim(0,21)
plt.xlim(0,80)
plt.xlabel('Energy (units: $\epsilon$)')
plt.ylabel('Ground state occupation, $n_0$')
plt.title('Manuscript Fig. 10 (upper)');
In [23]:
plt.figure()
plt.scatter(t_B,pop0_B[1:len(pop0_B)-1], label='bosons', color='blue')
plt.scatter(t_D,pop0_D[1:len(pop0_B)-1], label='distinguishable', color='red')
plt.xlim(0,2.5)
plt.ylim(0,21);
plt.xlabel('Temperature (units: $\epsilon/k$)')
plt.ylabel('Ground state occupation, $n_0$')
plt.title('Manuscript Fig. 10 (lower)')
plt.legend();
#plt.scatter(x,y)
In [26]:
plt.figure()
x1 = t_B
x2 = t_D
y = np.linspace(1,q_max-1, q_max-1)/n
plt.scatter(x1, y,color='blue', label='bosons')
plt.scatter(x2, y, color='red', label='distinguishable')
plt.xlim(0,2.5)
plt.ylim(0,5);
plt.xlabel('Temperature (units: $\epsilon/k$)')
plt.ylabel('Energy/particle (units: $\epsilon$)')
plt.title('Manuscript Fig. 11')
plt.legend(loc='upper left');
Example with more bosons: $N_B = 30$¶
In [139]:
n = 30
q_max = 90
pop0_B = np.zeros(q_max+1)
pop0_B[0] = n
s_B = np.zeros(q_max+1)
for q in range(1,q_max+1):
ex0_B = 0
omega_B = 0
for p in iterables.partitions(q,n):
npa = np.array(list(p.items())).T
#print(npa)
omega1 = mult_B(npa)
ex0_B += omega1*(n - np.sum(npa[1]))
omega_B += omega1
pop0_B[q] = ex0_B/omega_B
s_B[q] = np.log(omega_B)
print(q) # progress monitor
#print(ex0/omega, np.log(omega))
In [11]:
#pop0_B, s_B = np.loadtxt('temp').Tq
In [12]:
pop0_B
Out[12]:
In [13]:
t_B = 2/(s_B[2:] - s_B[:len(s_B)-2])
#t_D = 2/(s_D[2:] - s_D[:len(s_D)-2])
In [15]:
plt.figure()
plt.scatter(t_B,pop0_B[1:len(pop0_B)-1])
#plt.scatter(t_D,pop0_D[1:len(pop0_D)-1], color='red')
plt.xlim(0,2.5)
plt.ylim(0,26);
plt.xlabel('Temperature (units: $\epsilon/k$)')
plt.ylabel('Ground state occupation, $n_0$');
#plt.scatter(x,y)
Version information¶
version_information is from J.R. Johansson (jrjohansson at gmail.com); see Introduction to scientific computing with Python for more information and instructions for package installation.
In [2]:
%load_ext version_information
In [3]:
version_information numpy, scipy, matplotlib, sympy
Out[3]:
In [ ]: