Negative temperature and the old PHYS 211 entropy lab

In days of old in PHYS 211, before Ben's thermal physics unit, the entropy lab was the pegs and washers game. As a reminder, there were $N$ pegs around the circumference of circular board, which represented distinguishable particles with equally spaced energy levels. Indistinguishable units of energy were represented by washer, and the game started with the energy equally distributed, with one (or two) washers on each peg. Attached to the center of the board were two spinners, and an energy exchange was chosen by spinning them. One spinner pointed to the donor particle, and the other to the acceptor. Every 10 spins or so the distribution of energy (i.e., the number of particles with no washers ($n_0$), the number of particles with 1 washer ($n_1$), etc.) so that the evolution of the entropy could be monitored. The entropy in the old exercise is not the entropy used in the current PHYS 211 lab. In the current lab we're answering the question about which way heat flows between two systems; in the old lab we were investigating the distribution of energy within a closed system. In the old lab the multiplicity is

$$ W = \frac{N!}{n_0!\, n_1!\, n_2!\,\dots} $$

In this notebook I model a modified version of the pegs and washer game. The modification is that there are only a finite number of energy levels. In practice this means that once a peg has the maximum number of washers it isn't a valid acceptor.

In what follows I model a 5 level system, with energies 0, 1, 2, 3, and 4. I first investigate a system that starts with 2 washers per peg, and then investigate a system that starts with 4 washers per peg.

Comments on the python

• pegs 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 pegs to get the distribution
• the entropy 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

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 entropy(dist):
    '''Calculate entropy'''
    return np.log( special.factorial(n_pegs)/np.prod(special.factorial(dist)) )

Positive temperature example

Midpoint of energy scale is 2. In this example the energy per particle is less than 2, so the system doesn't exhibit population inversion.

n_pegs = 100    # Number of pegs
n_wpp = 1       # Initioal number of washers per peg
n_spins = 5000  # Number of energy exchanges


pegs = n_wpp*np.ones(n_pegs, dtype=int)  # Initialize array containing the number of washers on each peg
ent = np.array([entropy(np.bincount(pegs))])

# Sequence of energy exchanges and entropy calculations
i = 0

while i < n_spins:
    fr = stats.randint.rvs(0, n_pegs, size = 1)
    to = stats.randint.rvs(0, n_pegs, size = 1)
    if pegs[fr] > 0 and pegs[to] < n_levels - 1 and fr != to:
        pegs[fr] += -1
        pegs[to] += 1
        i += 1
        if i/10 == int(i/10):
            dist = np.bincount(pegs)
            ent = np.append(ent, entropy(dist))

dist = np.bincount(pegs)

plt.figure()
x = np.linspace(0,len(dist)-1,len(dist))
plt.bar(x, dist )
plt.xlabel('Number of washers')
plt.ylabel('Number of pegs');

plt.figure()
x = np.linspace(1, len(ent), len(ent))
y = ent
plt.plot(10*x, y)
plt.xlim(0,10*len(ent))
plt.xlabel('Exchange number')
plt.ylabel('Entropy')
plt.title('5 level system; intially 1 washer on each peg');

Negative temperature example

Midpoint of energy scale is 2. In this example the energy per particle is greater than 2, so the system does exhibit population inversion.

n_pegs = 100    # Number of pegs
n_wpp = 3       # Initioal number of washers per peg
n_spins = 5000  # Number of energy exchanges


pegs = n_wpp*np.ones(n_pegs, dtype=int)  # Initialize array containing the number of washers on each peg
ent = np.array([entropy(np.bincount(pegs))])

# Sequence of energy exchanges and entropy calculations
i = 0

while i < n_spins:
    fr = stats.randint.rvs(0, n_pegs, size = 1)
    to = stats.randint.rvs(0, n_pegs, size = 1)
    if pegs[fr] > 0 and pegs[to] < n_levels - 1 and fr != to:
        pegs[fr] += -1
        pegs[to] += 1
        i += 1
        if i/10 == int(i/10):
            dist = np.bincount(pegs)
            ent = np.append(ent, entropy(dist))

dist = np.bincount(pegs)

plt.figure()
x = np.linspace(0,len(dist)-1,len(dist))
plt.bar(x, dist )
plt.xlabel('Number of washers')
plt.ylabel('Number of pegs');

plt.figure()
x = np.linspace(1, len(ent), len(ent))
y = ent
plt.plot(10*x, y)
plt.xlabel('Exchange number')
plt.ylabel('Entropy')
plt.xlim(0,10*len(ent))
plt.title('5 level system; intially 3 washers on each peg');

%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 Jan 30 12:42:01 2023 EST