Einstein solid population distribution

Consider system of $N$ quantum harmonic oscillators with $M$ units of energy, i.e., $U = M\hbar \omega$. The muliplicity of the system is

$$ \Omega(N,M) = \frac{(M + N-1)!}{M! (N-1)!}, $$

and

$$ S(N,M) = k\ln \Omega(N,M). $$

Now consider counting the number of microstates with $N$ and $M$ and a specified number of oscillators with 0 units of energy, $n_0$, a specified number with 1 unit of energy, $n_1$, etc. It's not hard to show that this number is

$$ W(N,M; n_0, \dots, n_M) = \frac{N!}{n_0!\, n_1!\,\dots N_M!}. $$

The fraction of all equally probable microstates with this set of occupation numbers is

$$ P(N,M; n_0, \dots, n_M) = \frac{W(N,M;n_0, n_1\dots, n_M)}{\Omega(N,M))} $$

The most likely set of occupation numbers is thus the set that gives the largest value of $W$, or equivalently, $\ln(W)$.

The strategy used in this notebook is to start with some distribution of energy. (It doesn't really matter what the initial distribution is, but I start with the same number of energy units in each oscillator, i.e., a distribution with $W = 1$.) The energy is "shuffled" between oscillators many times while keeping track of the maximum value of $W$ that has been observed, along with the set of occupation numbers that resulted in $W_{\rm max}$.

The "shuffling" method used in this notebook is extremely simple:

• the oscillators are numbered 1 to $N$
• two random integers between 1 and $N$ are generated, the first random number corresponds to a "donor" oscillator, the second to an "acceptor" oscillator
• $n_{\rm donor}\rightarrow (n_{\rm donor} - 1)$ and $n_{\rm acceptor}\rightarrow (n_{\rm acceptor} + 1)$
• repeat

There are probably better ways to search for the largest value of $W$ and associated distribution, but this method is extremely easy to code, and makes the point that information about the occupation number distribution can be recovered (at least in principle) from the micrcanonical ensemble.

Comments on the python

• state is a numpy array in which each element gives the number of washers on a corresponding peg
• np.bincount is a very convenient way to bin state to get the distribution dist
• the value of $W$ ('w') is trivial to calculate by taking advantage of the ability to broadcast operations over the array dist, along with np.prod which calculates the product of the elements of an array

import numpy as np
from scipy import special
from scipy import stats

import matplotlib as mpl
import matplotlib.pyplot as plt

# Following is an Ipython magic command that puts figures in notebook.
%matplotlib notebook
        
# M.L. modification of matplotlib defaults
# Changes can also be put in matplotlibrc file, 
# or effected using mpl.rcParams[]
plt.style.use('classic')
plt.rc('figure', figsize = (6, 4.5)) # Reduces overall size of figures
plt.rc('axes', labelsize=16, titlesize=14)
plt.rc('xtick', labelsize=12)
plt.rc('ytick', labelsize=12)
plt.rc('figure', autolayout = True)  # Adjusts supblot params for new size

def w(dist):
    '''Calculate number of microstates for 'n_osc' oscillators with occupation numbers in array 'dist' '''
    return np.log( special.factorial(n_osc)/np.prod(special.factorial(dist)) )

n_osc = 100     # Number of oscillators
n_epo = 2       # Initial number of units of Energy Per Oscillator
m = n_osc*n_epo # Total energy M
n_exchanges = 10000  # Number of energy exchanges
w_max = 0


state = n_epo*np.ones(n_osc, dtype=int)  # Initialize array containing the number of energy units in each oscillator
w_history = np.array([w(np.bincount(state, minlength = m))])
np.bincount(state, minlength=m)

array([  0,   0, 100,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,   0,
         0,   0,   0,   0,   0])

state

array([2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
       2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])

Sequence of energy exchanges and calculation of $W$

i = 0

while i < n_exchanges:
    fr = stats.randint.rvs(0, n_osc, size = 1)  # Pick donor
    to = stats.randint.rvs(0, n_osc, size = 1)  # Pick acceptor
    if state[fr] > 0 and fr != to:
        state[fr] += -1
        state[to] += 1
        i += 1
    
        if i/10 == int(i/10):
            #print(pegs)
            dist = np.bincount(state, minlength = m)
            w_value = w(dist)
            if w_value > w_max:
                w_max = w_value
                dist_max = dist
            w_history = np.append(w_history, w_value)

Calculation of occupation number probablities for canonical ensemble

Using temperature determined for microcanonical ensemble

kt = 1/np.log(1 + 1/n_epo) # kT from microcanonical
es = np.linspace(1, m, m)  # energies
z = np.sum(np.exp(-es/kt)) # partition function
probabilities = np.exp(-es/kt)/z

plt.figure()
x = np.linspace(0,len(dist)-1, len(dist))
plt.bar(x, dist_max)
plt.scatter(x, probabilities*n_osc)
plt.xlim(-0.5,5*n_epo +0.5)
plt.xlabel('Number of units of energy')
plt.ylabel('Number of oscillator ($n_i$)');

Comment on graph

Agreement with Boltzman distribution is pretty good, even for these small values of $N$ and $M$! Value of $n_0$ is systematically higher than it would be in Boltzmann distribution. This difference decreases as the average energy per particle ($M/N$) increases (but simulations with larger values of $M/N$ require more exchanges.

plt.figure()
x = np.linspace(1, len(w_history), len(w_history))
y = w_history
plt.plot(10*x, y)
x2 = np.array([1,x[-1]])
y2 = np.array([w_max,w_max])
plt.plot(10*x2,y2, label='$W_{max}$')
plt.xlim(0,10*len(w_history))
plt.xlabel('Exchange number')
plt.ylabel('$\ln(W)$')
plt.legend(loc='lower right');

w_max

173.1597911822521

%load_ext version_information

version_information numpy, scipy,  matplotlib

Software	Version
Python	3.7.15 64bit [GCC 11.2.0]
IPython	7.31.1
OS	Linux 4.9.0 9 amd64 x86_64 with debian 9.13
numpy	1.21.5
scipy	1.7.3
matplotlib	3.5.3
Mon May 01 12:46:42 2023 EDT