Entropy from Counting: Einstein Solids

Derive thermodynamics from microstate counting using two Einstein solids.

Introduction

30 Second Summary

Have you ever wondered why a hot cup of coffee always cools down to room temperature and never spontaneously heats back up? The answer is not a mysterious "law" imposed on nature. It is pure mathematics and probability at an enormous scale.

In this project, you will build a computational simulation that proves the second law of thermodynamics is just statistics in disguise. Starting from the simplest possible model of a solid, you will count microstates, watch entropy emerge from combinatorics, and derive temperature from scratch.

What You'll Build

A Jupyter notebook that computes and visualizes the statistical mechanics of two interacting Einstein solids, following the pedagogical arc of Schroeder's Thermal Physics Chapter 2.

By the end of this project, you'll have:

  • A working multiplicity function that counts microstates for any Einstein solid, verified against known combinatorial results.
  • Plots showing the multiplicity peak sharpening with system size, proving that the second law of thermodynamics is a statistical certainty.
  • Curves of entropy and temperature derived purely from counting, with a heat capacity plot showing both quantum freezeout and the classical Dulong-Petit limit.
  • Secret Mission: Simulate random energy exchanges between oscillators and watch the Boltzmann distribution emerge from pure randomness.

Are there any prerequisites?

You will need basic Python knowledge (loops, functions, arrays) and familiarity with factorials and combinatorics. Some calculus (derivatives) is helpful for the temperature and heat capacity sections.

No prior thermodynamics or statistical mechanics knowledge is required. The project builds everything from scratch using counting arguments.

Set Up Your Environment

Before you can count microstates and watch thermodynamics emerge, you need a computational workspace. JupyterLab is an interactive notebook environment where you can mix code, plots, and explanations in a single document.

In this step, you will install the scientific Python stack and create the notebook that will hold your entire computational exploration.

In this step, get ready to:
  • Install Python 3.11+, JupyterLab, NumPy, and Matplotlib.
  • Launch JupyterLab and create your project notebook.
  • Add a markdown header and import the libraries you will use throughout the project.
Install Python and JupyterLab

This project requires Python 3.11 or higher. NumPy 2.4+ depends on it, and you will use features like math.comb that require at least Python 3.8.

  • Press Cmd+Space (macOS) or the Windows key (Windows) to open your search bar.
  • Type Terminal (macOS) or PowerShell (Windows) and press Enter to open it.
  • Check your Python version by running this command:
python3 --version

✔️ I see version 3.11 or higher

Your Python version is ready. Move on to the next instruction below.

ⓧ I see an older version

macOS ships with an older system Python that may shadow newer installs. You need to install Python 3.11+ side by side.

  • Download the latest Python 3.11+ installer from python.org/downloads.
  • Run the installer and follow the prompts.
  • After installation, close and reopen your terminal, then run python3 --version again to confirm it shows 3.11 or higher.

macOS path priority

If python3 --version still shows the old version after installing, the system Python is shadowing the new one. Try running /usr/local/bin/python3 --version or use the full path from the installer's Applications/Python 3.x/ folder.

ⓧ Command not found

Python is not installed yet. Install it now.

  • Go to python.org/downloads and download the latest Python 3.11+ installer for your operating system.
  • Run the installer and follow the prompts to completion.
  • Close and reopen your terminal.
  • Run python3 --version to confirm the installation shows 3.11 or higher.
  • Move to your Desktop by running this command:
cd ~/Desktop
  • Create your project folder and move into it by running these commands:
mkdir einstein-solids
cd einstein-solids
  • Install JupyterLab and the scientific Python packages by running:
pip install jupyterlab numpy matplotlib scipy

What does this command do?

  • jupyterlab gives you the interactive notebook environment.
  • numpy provides fast array operations for computing multiplicities across all macrostates.
  • matplotlib handles all the plotting (multiplicity peaks, entropy curves, heat capacity).
  • scipy provides gammaln for computing log-factorials without integer overflow on large systems.

Seeing 'pip: command not found'?

  • Try pip3 install jupyterlab numpy matplotlib scipy instead.
  • If that also fails, run python3 -m pip install jupyterlab numpy matplotlib scipy to invoke pip through your Python 3 installation directly.

Help me fix my pip installation

  • Verify the installations by running these commands:
jupyter lab --version
python3 -c "import numpy; print(numpy.__version__)"
python3 -c "import matplotlib; print(matplotlib.__version__)"

You should see version numbers printed for each command with no errors.

  • Launch JupyterLab by running:
jupyter lab

JupyterLab opens in your browser, showing a Launcher tab with options to create new files.

  • In the Launcher, click Python 3 under the Notebook section to create a new notebook.
  • Right-click the Untitled.ipynb tab at the top and select Rename Notebook.
  • Type einstein_solids.ipynb and press Enter.

What is a Jupyter notebook?

A Jupyter notebook is a document made of cells. Each cell is either code (which you can run) or markdown (formatted text for explanations). You will use both throughout this project to combine computation with narrative.

Set up the notebook with imports and initial markdown

Every good computational physics notebook starts with a title and the libraries it depends on. You will add a markdown header describing the project, then an import cell that loads all the tools you need.

  • The first cell in your notebook is a code cell by default. Change it to a markdown cell by pressing Esc to enter command mode, then pressing M on your keyboard.
  • Type the following markdown into the cell:
# Entropy from Counting: Two Einstein Solids in Thermal Contact

Following Schroeder's *Thermal Physics* Chapter 2 approach:
start by counting microstates, then watch thermodynamics emerge.

What does this markdown do?

The # creates a top-level heading. The *Thermal Physics* wraps the book title in italics. When you run the cell, the raw markdown renders into formatted text.

  • Run the markdown cell by pressing Shift+Enter.

The cell renders into a formatted heading and paragraph. A new empty code cell appears below it.

  • In the new code cell (Cell 2), type the following imports:
import numpy as np
import matplotlib.pyplot as plt
from math import comb, log
from scipy.special import gammaln

What does this code do?

  • numpy as np gives you fast array operations for computing multiplicities across hundreds of macrostates at once.
  • matplotlib.pyplot as plt provides all the plotting functions for visualizing distributions and curves.
  • comb from math computes exact binomial coefficients (the core of the multiplicity formula).
  • log from math computes natural logarithms for entropy calculations.
  • gammaln from scipy.special computes log-factorials without overflow, essential for large systems where regular factorials exceed integer limits.
  • Run the import cell by pressing Shift+Enter.
  • Confirm the cell runs with no errors. You should see a small number (like [1] or [2]) appear to the left of the cell, indicating successful execution.

Seeing a ModuleNotFoundError?

  • Make sure you launched jupyter lab from the same terminal session where you ran pip install.
  • If you have multiple Python versions, the notebook kernel might be using a different one. Go to Kernel > Change Kernel and select the Python 3.11+ kernel.

Help me fix my import errors in JupyterLab

✔️ Awesome, I've got everything!

Great. Make sure you have saved your notebook by pressing Cmd+S (macOS) or Ctrl+S (Windows).

ⓧ I'd like to double check the full code

Here is what your notebook should contain at this point.

Cell 1 (Markdown):

# Entropy from Counting: Two Einstein Solids in Thermal Contact

Following Schroeder's *Thermal Physics* Chapter 2 approach:
start by counting microstates, then watch thermodynamics emerge.

Cell 2 (Code):

import numpy as np
import matplotlib.pyplot as plt
from math import comb, log
from scipy.special import gammaln

Your computational workspace is ready. Next up, you will write the multiplicity function and start counting microstates for the Einstein solid.

Compute Multiplicity of a Single Einstein Solid

Your environment is ready and your notebook has its imports loaded.

Now you get to the physics. The Einstein solid is the simplest model where thermodynamics emerges from pure counting. It consists of N quantum harmonic oscillators sharing q indistinguishable energy quanta. Your first job is to compute how many ways those quanta can be distributed.

In this step, get ready to:
  • Implement the multiplicity function using the stars-and-bars formula.
  • Verify your function against known results for small systems.
  • Plot multiplicity growth to see factorial-scale explosion.
Implement the multiplicity function

The multiplicity of an Einstein solid counts the total number of microstates (ways to distribute q energy quanta among N oscillators). The formula comes from combinatorics, specifically the "stars and bars" technique.

  • Create a new Markdown cell below your import cell. Add the following section header and formula explanation:
## Step 2: Multiplicity of a Single Einstein Solid

An Einstein solid has N quantum harmonic oscillators sharing q indistinguishable
energy quanta. The multiplicity (number of microstates) is the "stars and bars"
combinatorial formula:

$$\Omega(N, q) = \binom{q + N - 1}{q}$$

What is the stars and bars formula?

Imagine q identical balls (energy quanta) separated by N-1 dividers (the walls between oscillators). The number of ways to arrange these is C(q + N - 1, q). This is a standard combinatorial identity.

  • Create a new code cell below the markdown cell. Define the multiplicity function and verify it against known values by entering this code:
def multiplicity(N, q):
    """Number of microstates for N oscillators sharing q energy quanta."""
    return comb(q + N - 1, q)

# Verify against Schroeder's examples
print(f"Omega(4, 10) = {multiplicity(4, 10)}")  # Should be 286
print(f"Omega(4, 2) = {multiplicity(4, 2)}")    # Should be 10
print(f"Omega(3, 3) = {multiplicity(3, 3)}")    # Should be 10

What does this code do?

The multiplicity function uses Python's math.comb to compute the exact combinatorial result C(q + N - 1, q). This returns an exact Python integer with arbitrary precision, which is important for verifying small systems by hand.

The verification values come from Schroeder's textbook. For example, 4 oscillators sharing 10 quanta gives C(13, 10) = 286 microstates.

  • Run the cell by pressing Shift+Enter.
  • You should see three printed lines: Omega(4, 10) = 286, Omega(4, 2) = 10, and Omega(3, 3) = 10.

Getting a NameError for comb?

Make sure you ran Cell 2 (your import cell) first. Jupyter cells must be executed in order for variables and functions to be available.

If you restarted the kernel at any point, re-run Cell 2 before running Cell 4.

Help me debug my multiplicity function

Explore small systems by hand

Before scaling up to large systems, it helps to build intuition by examining several small configurations. This lets you spot-check results against what you could count on paper.

  • Create a new code cell below Cell 4. Compute multiplicity for a range of small systems by entering this code:
# Explore small systems
for N in [3, 4, 10]:
    for q in [2, 5, 10]:
        print(f"N={N}, q={q}: Omega = {multiplicity(N, q)}")

What does this code do?

This nested loop iterates over combinations of oscillator counts (N) and energy quanta (q). For each combination, it prints the multiplicity so you can see how the number of microstates grows with both N and q.

Notice how even for N=10 and q=10, the multiplicity reaches 92,378. The growth is combinatorial, not linear.

  • Run the cell with Shift+Enter.
  • You should see 9 printed lines showing multiplicity values. Verify that N=3, q=3 gives Omega = 10, which you can confirm by hand as C(5,3) = 10.

Seeing unexpected values?

Double check that your multiplicity function uses q + N - 1 as the first argument to comb, not N + q - 1. Both are mathematically equivalent, but a typo like N - 1 in the second argument would give wrong results.

Help me fix my multiplicity values

Plot multiplicity growth

Numbers in a table only tell part of the story. A plot reveals just how explosively multiplicity grows. You will need a logarithmic y-axis because the values span many orders of magnitude.

  • Create a new code cell below Cell 5. Plot multiplicity vs. energy quanta for a fixed number of oscillators by entering this code:
# Plot multiplicity growth for fixed N
N = 50
q_values = np.arange(0, 101)
omega_values = [multiplicity(N, q) for q in q_values]

plt.figure(figsize=(8, 5))
plt.semilogy(q_values, omega_values)
plt.xlabel("Energy quanta q")
plt.ylabel("Multiplicity (log scale)")
plt.title(f"Multiplicity of an Einstein Solid (N={N} oscillators)")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

What does this code do?

This computes multiplicity for N=50 oscillators as the energy q ranges from 0 to 100. The plt.semilogy function creates a plot with a logarithmic y-axis. Without the log scale, you would only see a vertical spike at high q values because the growth is so extreme.

Even on a log scale, the curve shoots upward. This factorial-scale growth is what makes large systems so predictable. The system overwhelmingly prefers the macrostate with the most microstates.

  • Run the cell with Shift+Enter.
  • You should see a plot with an upward-curving line on a log-scale y-axis, showing multiplicity rising from 1 (at q=0) to astronomically large values (above 10^30 at q=100).
  • Save your notebook with Cmd+S (macOS) or Ctrl+S (Windows).

Plot not showing up?

Make sure you ran your import cell (Cell 2) that imports matplotlib.pyplot. If the plot opens in a separate window instead of inline, add %matplotlib inline at the top of your import cell and re-run it.

Help me get my matplotlib plot to display

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook before moving on.

ⓧ I'd like to double check the full code

Here is the complete state of your notebook through the end of this step. Each section below represents one cell.

# Entropy from Counting: Two Einstein Solids in Thermal Contact

Following Schroeder's *Thermal Physics* Chapter 2 approach:
start by counting microstates, then watch thermodynamics emerge.
import numpy as np
import matplotlib.pyplot as plt
from math import comb, log
from scipy.special import gammaln
## Step 2: Multiplicity of a Single Einstein Solid

An Einstein solid has N quantum harmonic oscillators sharing q indistinguishable
energy quanta. The multiplicity (number of microstates) is the "stars and bars"
combinatorial formula:

$$\Omega(N, q) = \binom{q + N - 1}{q}$$
def multiplicity(N, q):
    """Number of microstates for N oscillators sharing q energy quanta."""
    return comb(q + N - 1, q)

# Verify against Schroeder's examples
print(f"Omega(4, 10) = {multiplicity(4, 10)}")  # Should be 286
print(f"Omega(4, 2) = {multiplicity(4, 2)}")    # Should be 10
print(f"Omega(3, 3) = {multiplicity(3, 3)}")    # Should be 10
# Explore small systems
for N in [3, 4, 10]:
    for q in [2, 5, 10]:
        print(f"N={N}, q={q}: Omega = {multiplicity(N, q)}")
# Plot multiplicity growth for fixed N
N = 50
q_values = np.arange(0, 101)
omega_values = [multiplicity(N, q) for q in q_values]

plt.figure(figsize=(8, 5))
plt.semilogy(q_values, omega_values)
plt.xlabel("Energy quanta q")
plt.ylabel("Multiplicity (log scale)")
plt.title(f"Multiplicity of an Einstein Solid (N={N} oscillators)")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

You now have a verified multiplicity function and visual proof of how explosively it grows. Next up, you will put two Einstein solids in thermal contact and watch where thermodynamics comes from.

Simulate Two Solids in Thermal Contact

You've verified that the multiplicity of a single Einstein solid grows astronomically with energy. A single solid sitting alone isn't very interesting physically, though.

The real question is: when two solids can exchange energy, which energy split is the system most likely to be found in? In this step, you'll put two solids in thermal contact and compute the combined multiplicity for every possible way to divide energy between them.

In this step, get ready to:
  • Compute the combined multiplicity across all macrostates of two interacting solids.
  • Plot the multiplicity peak for small and large systems.
  • Quantify how fluctuations vanish as system size grows.
Compute the combined multiplicity for all macrostates

Two Einstein solids A and B share a fixed total energy q_total. For each possible split (giving q_A quanta to solid A and the rest to solid B), the combined multiplicity is the product of the individual multiplicities.

What is a macrostate?

A macrostate describes the large-scale configuration of the system. Here, it's simply the energy split: how many quanta solid A has vs. solid B.

Each macrostate contains many microstates (specific arrangements of quanta among oscillators). The macrostate with the most microstates is overwhelmingly the most likely to be observed.

  • Create a new Markdown cell below Cell 6 and add this section header:
## Step 3: Two Solids in Thermal Contact

Two solids A and B share a fixed total energy $q_{total}$.
The combined multiplicity for a given split is:

$$\Omega_{total}(q_A) = \Omega_A(q_A) \times \Omega_B(q_{total} - q_A)$$
  • Create a new code cell below the markdown cell and define the combined multiplicity function:
def combined_multiplicity(N_A, N_B, q_total):
    """Compute combined multiplicity for each energy split."""
    q_A_values = np.arange(0, q_total + 1)
    omega_combined = np.array([
        multiplicity(N_A, q_A) * multiplicity(N_B, q_total - q_A)
        for q_A in q_A_values
    ])
    return q_A_values, omega_combined

What does this code do?

For every possible value of q_A (from 0 to q_total), the function computes the product of the multiplicity of solid A (with q_A quanta) and solid B (with the remaining quanta).

This product is the total number of microstates accessible to the composite system for that particular macrostate. The function returns both the array of q_A values and the corresponding combined multiplicities.

  • Run the cell.

You won't see any output yet. The next cell will call this function and plot the results.

Plot the multiplicity peak for small and large systems

With the function defined, you can now visualize which energy splits are most likely. You'll start with a small system where the peak is broad, then scale up to watch it sharpen dramatically.

  • Create a new code cell and plot the combined multiplicity for a small system (N_A = N_B = 10, q_total = 20):
# Small system: N_A = N_B = 10, q_total = 20
q_A, omega = combined_multiplicity(10, 10, 20)
probability = omega / omega.sum()

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].plot(q_A, omega)
axes[0].set_xlabel("Energy in solid A (q_A)")
axes[0].set_ylabel("Combined multiplicity")
axes[0].set_title("N_A = N_B = 10, q_total = 20")
axes[0].grid(True, alpha=0.3)

axes[1].plot(q_A, probability)
axes[1].set_xlabel("Energy in solid A (q_A)")
axes[1].set_ylabel("Probability P(q_A)")
axes[1].set_title("Probability distribution")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Most probable q_A: {q_A[np.argmax(omega)]}")
print(f"Peak probability: {probability.max():.4f}")

What does this code do?

  • It calls combined_multiplicity for two 10-oscillator solids sharing 20 quanta.
  • It normalizes the combined multiplicity to get a probability distribution (dividing by the total number of microstates across all macrostates).
  • The left plot shows raw multiplicity vs. energy split. The right plot shows the same data as a probability.
  • Run the cell.

You should see a broad peak centered at q_A = 10 (equal energy sharing for equal-sized solids). The peak probability is only about 0.07, meaning the system has a fair chance of being found far from equilibrium.

Not seeing a peak at q_A = 10?

  • Make sure Cell 8 (the combined_multiplicity function) ran without errors before running Cell 9.
  • Confirm that Cell 4 (the multiplicity function) also ran successfully. Cell 9 depends on both functions being defined.

Why isn't my combined multiplicity plot showing a peak?

Now you'll see what happens when the system gets larger. This is where the physics gets dramatic.

  • Create a new code cell and plot the probability distribution for increasing system sizes:
# Large system: observe the sharpening
fig, ax = plt.subplots(figsize=(8, 5))

for N in [10, 50, 100, 200]:
    q_total = 2 * N  # Equal average energy per oscillator
    q_A_vals, omega_vals = combined_multiplicity(N, N, q_total)
    prob = omega_vals / omega_vals.sum()
    # Normalize x-axis to fraction of total energy
    ax.plot(q_A_vals / q_total, prob, label=f"N = {N}")

ax.set_xlabel("Fraction of energy in solid A (q_A / q_total)")
ax.set_ylabel("Probability")
ax.set_title("Sharpening of the multiplicity peak with system size")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

What does this code do?

  • It loops over four system sizes (N = 10, 50, 100, 200), keeping the average energy per oscillator fixed at 2 quanta.
  • The x-axis is normalized to the fraction of total energy in solid A, so all curves are comparable on the same scale.
  • As N increases, the probability distribution gets sharper and taller. For N = 200, the system is virtually certain to be found within a few percent of equal energy sharing.
  • Run the cell.

You should see four curves overlaid, each sharper than the last. The N = 200 curve is so narrow it looks like a spike compared to the broad N = 10 curve.

Why does this matter?

This sharpening IS the second law of thermodynamics in action. For a real solid with ~10^23 oscillators, the peak is so fantastically sharp that the system is never observed away from the most probable macrostate.

The second law isn't a separate axiom. It's a statistical certainty that emerges from counting.

Cell 10 taking a long time?

For N = 200, the loop computes 401 multiplicity values for each solid. This should finish in under a second on most machines. If it hangs, check that you haven't accidentally passed a much larger N value.

Why is my Einstein solid computation running slowly?

Quantify the sharpness

The plot shows the peak getting sharper, but you can be more precise. By computing the standard deviation of the probability distribution at each system size, you can show that relative fluctuations decrease as 1/sqrt(N).

  • Create a new code cell and compute the relative fluctuation as a function of system size:
# Quantify: standard deviation vs system size
sizes = [10, 20, 50, 100, 200, 500]
std_devs = []
relative_flucts = []

for N in sizes:
    q_total = 2 * N
    q_A_vals, omega_vals = combined_multiplicity(N, N, q_total)
    prob = omega_vals / omega_vals.sum()
    mean = np.sum(q_A_vals * prob)
    variance = np.sum((q_A_vals - mean)**2 * prob)
    std = np.sqrt(variance)
    std_devs.append(std)
    relative_flucts.append(std / mean)

plt.figure(figsize=(8, 5))
plt.loglog(sizes, relative_flucts, 'o-')
plt.loglog(sizes, 1.0 / np.sqrt(sizes), '--', label="1/sqrt(N) scaling")
plt.xlabel("System size N")
plt.ylabel("Relative fluctuation (sigma / mean)")
plt.title("Fluctuations vanish for large systems")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

What does this code do?

  • For each system size, it computes the mean and standard deviation of the probability distribution over macrostates.
  • The relative fluctuation (sigma / mean) tells you what fraction of the total energy the system typically deviates from equilibrium.
  • The log-log plot confirms the 1/sqrt(N) scaling. The data points should fall right on the dashed reference line.
  • Run the cell.

You should see the data points falling along the dashed 1/sqrt(N) line. For N = 500, the relative fluctuation is about 3%, meaning the system stays within 3% of equal energy sharing.

What does 1/sqrt(N) scaling mean physically?

For a real solid with N ~ 10^23 oscillators, the relative fluctuation is about 10^(-12). The most probable macrostate IS effectively the only macrostate. This is why you never see energy spontaneously flow from cold to hot.

Data points not following the 1/sqrt(N) line?

  • Make sure your sizes list includes a good range. The scaling is approximate for very small N (like 10) and becomes more exact for larger values.
  • Check that you're dividing standard deviation by the mean (not by N or q_total). The relative fluctuation is sigma / mean.

Why aren't my fluctuation data points following the expected 1/sqrt(N) scaling?

✔️ Awesome, I've got everything!

Great. Make sure you've saved your notebook before moving on.

ⓧ I'd like to double check the full code

Here is the complete notebook up to this point. Compare against your cells:

# Cell 1 (Markdown)
# # Entropy from Counting: Two Einstein Solids in Thermal Contact
# 
# Following Schroeder's *Thermal Physics* Chapter 2 approach:
# start by counting microstates, then watch thermodynamics emerge.

# Cell 2
import numpy as np
import matplotlib.pyplot as plt
from math import comb, log
from scipy.special import gammaln

# Cell 3 (Markdown)
# ## Step 2: Multiplicity of a Single Einstein Solid
#
# An Einstein solid has N quantum harmonic oscillators sharing q indistinguishable
# energy quanta. The multiplicity (number of microstates) is the "stars and bars"
# combinatorial formula:
#
# $$\Omega(N, q) = \binom{q + N - 1}{q}$$

# Cell 4
def multiplicity(N, q):
    """Number of microstates for N oscillators sharing q energy quanta."""
    return comb(q + N - 1, q)

# Verify against Schroeder's examples
print(f"Omega(4, 10) = {multiplicity(4, 10)}")  # Should be 286
print(f"Omega(4, 2) = {multiplicity(4, 2)}")    # Should be 10
print(f"Omega(3, 3) = {multiplicity(3, 3)}")    # Should be 10

# Cell 5
# Explore small systems
for N in [3, 4, 10]:
    for q in [2, 5, 10]:
        print(f"N={N}, q={q}: Omega = {multiplicity(N, q)}")

# Cell 6
# Plot multiplicity growth for fixed N
N = 50
q_values = np.arange(0, 101)
omega_values = [multiplicity(N, q) for q in q_values]

plt.figure(figsize=(8, 5))
plt.semilogy(q_values, omega_values)
plt.xlabel("Energy quanta q")
plt.ylabel("Multiplicity (log scale)")
plt.title(f"Multiplicity of an Einstein Solid (N={N} oscillators)")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Cell 7 (Markdown)
# ## Step 3: Two Solids in Thermal Contact
#
# Two solids A and B share a fixed total energy $q_{total}$.
# The combined multiplicity for a given split is:
#
# $$\Omega_{total}(q_A) = \Omega_A(q_A) \times \Omega_B(q_{total} - q_A)$$

# Cell 8
def combined_multiplicity(N_A, N_B, q_total):
    """Compute combined multiplicity for each energy split."""
    q_A_values = np.arange(0, q_total + 1)
    omega_combined = np.array([
        multiplicity(N_A, q_A) * multiplicity(N_B, q_total - q_A)
        for q_A in q_A_values
    ])
    return q_A_values, omega_combined

# Cell 9
# Small system: N_A = N_B = 10, q_total = 20
q_A, omega = combined_multiplicity(10, 10, 20)
probability = omega / omega.sum()

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].plot(q_A, omega)
axes[0].set_xlabel("Energy in solid A (q_A)")
axes[0].set_ylabel("Combined multiplicity")
axes[0].set_title("N_A = N_B = 10, q_total = 20")
axes[0].grid(True, alpha=0.3)

axes[1].plot(q_A, probability)
axes[1].set_xlabel("Energy in solid A (q_A)")
axes[1].set_ylabel("Probability P(q_A)")
axes[1].set_title("Probability distribution")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Most probable q_A: {q_A[np.argmax(omega)]}")
print(f"Peak probability: {probability.max():.4f}")

# Cell 10
# Large system: observe the sharpening
fig, ax = plt.subplots(figsize=(8, 5))

for N in [10, 50, 100, 200]:
    q_total = 2 * N  # Equal average energy per oscillator
    q_A_vals, omega_vals = combined_multiplicity(N, N, q_total)
    prob = omega_vals / omega_vals.sum()
    # Normalize x-axis to fraction of total energy
    ax.plot(q_A_vals / q_total, prob, label=f"N = {N}")

ax.set_xlabel("Fraction of energy in solid A (q_A / q_total)")
ax.set_ylabel("Probability")
ax.set_title("Sharpening of the multiplicity peak with system size")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Cell 11
# Quantify: standard deviation vs system size
sizes = [10, 20, 50, 100, 200, 500]
std_devs = []
relative_flucts = []

for N in sizes:
    q_total = 2 * N
    q_A_vals, omega_vals = combined_multiplicity(N, N, q_total)
    prob = omega_vals / omega_vals.sum()
    mean = np.sum(q_A_vals * prob)
    variance = np.sum((q_A_vals - mean)**2 * prob)
    std = np.sqrt(variance)
    std_devs.append(std)
    relative_flucts.append(std / mean)

plt.figure(figsize=(8, 5))
plt.loglog(sizes, relative_flucts, 'o-')
plt.loglog(sizes, 1.0 / np.sqrt(sizes), '--', label="1/sqrt(N) scaling")
plt.xlabel("System size N")
plt.ylabel("Relative fluctuation (sigma / mean)")
plt.title("Fluctuations vanish for large systems")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

You've now seen the second law emerge from pure counting. The most probable macrostate becomes overwhelmingly dominant as systems grow. Next, you'll turn this multiplicity into entropy and temperature, showing that thermal equilibrium is just entropy maximization in disguise.

Derive Entropy and Temperature from Counting

In the previous step, you proved that large systems are overwhelmingly likely to be found near one specific energy split. Now the question becomes: what physical quantity tells us which split?

The answer is entropy. In this step, you will define entropy as the logarithm of multiplicity, then derive temperature as the slope of that entropy curve. The punchline: the second law of thermodynamics is not a separate axiom. It falls out of pure combinatorics.

In this step, get ready to:
  • Define entropy using the log of multiplicity and handle large-system overflow with gammaln.
  • Derive temperature from the numerical derivative of the entropy curve.
  • Verify that entropy maximization and equal temperatures are the same physical condition.
Define entropy and compute it

Entropy is defined as S = k * ln(Omega), where Omega is the multiplicity. Working in natural units (k = 1) simplifies this to S = ln(Omega).

For large systems, computing multiplicity(N, q) directly produces astronomically large integers that overflow even Python's arbitrary-precision arithmetic when you try to take their log. The solution is to compute ln(Omega) directly using gammaln from SciPy.

  • Add a new markdown cell in your notebook to introduce this section:
# Cell 12 (Markdown)
# ## Step 4: Entropy and Temperature from Counting
#
# Entropy: $S = k \ln \Omega$ (using $k = 1$ natural units)
#
# Temperature: $\frac{1}{T} = \frac{\partial S}{\partial U}$ where $U = q \epsilon$

What does this markdown say?

This cell documents the two key formulas you will implement. Entropy is the log of multiplicity, and temperature is defined through the partial derivative of entropy with respect to energy.

  • In a new code cell below the markdown, define the log_multiplicity function and verify it against the exact result:
def log_multiplicity(N, q):
    """Compute ln(Omega) using gammaln to avoid overflow for large systems."""
    return gammaln(q + N) - gammaln(q + 1) - gammaln(N)

# Verify against exact for small systems
print(f"Exact: ln(Omega(4,10)) = {log(multiplicity(4, 10)):.6f}")
print(f"Approx: ln(Omega(4,10)) = {log_multiplicity(4, 10):.6f}")

What does this code do?

The gammaln function computes the logarithm of the gamma function. Since Gamma(n+1) = n!, this lets you compute ln(n!) without ever forming the huge integer n! itself.

The formula gammaln(q + N) - gammaln(q + 1) - gammaln(N) computes ln(C(q+N-1, q)) using the log-subtraction identity ln(a/b) = ln(a) - ln(b). The verification confirms both methods agree for small systems.

  • Run the cell.
  • Both printed values should match to six decimal places, confirming that log_multiplicity reproduces the exact result.

Seeing a NameError for gammaln?

Make sure Cell 2 (your imports cell) includes from scipy.special import gammaln. If you added it after running Cell 2 the first time, re-run Cell 2 before running this cell.

Help me fix my gammaln import error.

Now you have a function that computes entropy (as ln(Omega)) for systems of any size without overflow. Time to use it on the two-solid system.

  • In a new code cell, compute the entropy of both solids across all energy splits and plot the result:
# Entropy of two solids as a function of q_A
N_A, N_B, q_total = 100, 100, 200

q_A_vals = np.arange(0, q_total + 1)
S_A = np.array([log_multiplicity(N_A, q) for q in q_A_vals])
S_B = np.array([log_multiplicity(N_B, q_total - q) for q in q_A_vals])
S_total = S_A + S_B

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(q_A_vals, S_A, label="S_A (solid A)")
ax.plot(q_A_vals, S_B, label="S_B (solid B)")
ax.plot(q_A_vals, S_total, label="S_total = S_A + S_B", linewidth=2)
ax.axvline(q_A_vals[np.argmax(S_total)], color='red', linestyle='--',
           label=f"Equilibrium: q_A = {q_A_vals[np.argmax(S_total)]}")
ax.set_xlabel("Energy in solid A (q_A)")
ax.set_ylabel("Entropy S/k")
ax.set_title("Entropy maximization determines equilibrium")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

What does this code do?

For each possible energy split (q_A from 0 to 200), the code computes the entropy of solid A, the entropy of solid B (which gets the remaining energy), and the total entropy S_A + S_B.

Since entropy is a logarithm, the total entropy of two independent systems is the sum of individual entropies (just as ln(A * B) = ln(A) + ln(B)). The red dashed line marks where S_total is maximized.

  • Run the cell. You should see three curves: S_A rising, S_B falling, and S_total with a clear peak at q_A = 100 (equal energy sharing for equal-sized solids).

Plot is blank or S_total has no visible peak?

Check that you defined log_multiplicity in a previous cell and ran it successfully. If the function is not defined, this cell will raise a NameError.

Also verify that q_total is set to 200 (not 20). A very small system will have a broad, flat peak that is hard to see.

Help me debug my entropy plot.

Derive temperature from the entropy curve

Temperature in statistical mechanics is defined by 1/T = dS/dU. In natural units with epsilon = 1, this simplifies to 1/T = dS/dq.

If you compute the numerical derivative of S_A and S_B with respect to q_A, the point where the two slopes are equal is where the two solids have the same temperature. That is thermal equilibrium.

  • In a new code cell, compute the numerical derivatives and plot the inverse temperatures:
# Derive temperature: 1/T = dS/dq (with epsilon = 1)
dS_A = np.diff(S_A)  # Numerical derivative
dS_B = -np.diff(S_B)  # Note: as q_A increases, q_B decreases

q_midpoints = q_A_vals[:-1] + 0.5

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(q_midpoints, dS_A, label="dS_A/dq_A = 1/T_A")
ax.plot(q_midpoints, dS_B, label="-dS_B/dq_A = 1/T_B")
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_xlabel("q_A")
ax.set_ylabel("1/T (inverse temperature)")
ax.set_title("Equilibrium: where temperatures are equal (curves cross)")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Find crossing point
crossing_idx = np.argmin(np.abs(dS_A - dS_B))
print(f"Temperatures equal at q_A ~ {q_midpoints[crossing_idx]:.1f}")
print(f"Maximum of S_total at q_A = {q_A_vals[np.argmax(S_total)]}")

What does this code do?

The np.diff(S_A) call computes the finite difference dS_A/dq_A, which approximates the derivative. This gives 1/T_A at each energy value.

For solid B, the sign is negated because as q_A increases, q_B decreases. So -np.diff(S_B) gives the inverse temperature of solid B expressed in terms of q_A.

Where the two curves cross, both solids share the same temperature. The code finds this crossing point and prints it alongside the entropy-maximum q_A for comparison.

  • Run the cell. You should see two curves (1/T_A decreasing, 1/T_B increasing) crossing near q_A = 100.
  • Check the printed output. Both values should show q_A = 100, confirming that equal temperatures and maximum entropy point to the same energy split.

Curves don't cross or printed values don't match?

If the curves never cross, check the sign on dS_B. It must be -np.diff(S_B) (negated) because solid B loses energy when solid A gains it.

If the printed values differ by more than 1, make sure that S_A, S_B, and S_total were computed in the previous cell with the same q_A_vals array. Re-run both cells in order.

Help me understand why my temperature curves aren't crossing.

Verify the equilibrium condition

Look at the two printed values from the cell you just ran. They should both read q_A = 100.

Why does this matter?

You just proved two statements are the same thing. "The system maximizes total entropy" (from counting microstates) and "the two subsystems reach the same temperature" (from the thermodynamic definition of T) both point to the identical energy split.

This is the second law of thermodynamics derived from pure combinatorics. Energy flows from hot to cold because that direction has overwhelmingly more microstates available. No new axiom required.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook.

ⓧ I'd like to double check the full code

Here is the full notebook up to this point. Compare your cells against this reference:

# Cell 1 (Markdown)
# # Entropy from Counting: Two Einstein Solids in Thermal Contact
# 
# Following Schroeder's *Thermal Physics* Chapter 2 approach:
# start by counting microstates, then watch thermodynamics emerge.

# Cell 2
import numpy as np
import matplotlib.pyplot as plt
from math import comb, log
from scipy.special import gammaln

# Cell 3 (Markdown)
# ## Step 2: Multiplicity of a Single Einstein Solid
#
# An Einstein solid has N quantum harmonic oscillators sharing q indistinguishable
# energy quanta. The multiplicity (number of microstates) is the "stars and bars"
# combinatorial formula:
#
# $$\Omega(N, q) = \binom{q + N - 1}{q}$$

# Cell 4
def multiplicity(N, q):
    """Number of microstates for N oscillators sharing q energy quanta."""
    return comb(q + N - 1, q)

# Verify against Schroeder's examples
print(f"Omega(4, 10) = {multiplicity(4, 10)}")  # Should be 286
print(f"Omega(4, 2) = {multiplicity(4, 2)}")    # Should be 10
print(f"Omega(3, 3) = {multiplicity(3, 3)}")    # Should be 10

# Cell 5
# Explore small systems
for N in [3, 4, 10]:
    for q in [2, 5, 10]:
        print(f"N={N}, q={q}: Omega = {multiplicity(N, q)}")

# Cell 6
# Plot multiplicity growth for fixed N
N = 50
q_values = np.arange(0, 101)
omega_values = [multiplicity(N, q) for q in q_values]

plt.figure(figsize=(8, 5))
plt.semilogy(q_values, omega_values)
plt.xlabel("Energy quanta q")
plt.ylabel("Multiplicity (log scale)")
plt.title(f"Multiplicity of an Einstein Solid (N={N} oscillators)")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Cell 7 (Markdown)
# ## Step 3: Two Solids in Thermal Contact
#
# Two solids A and B share a fixed total energy $q_{total}$.
# The combined multiplicity for a given split is:
#
# $$\Omega_{total}(q_A) = \Omega_A(q_A) \times \Omega_B(q_{total} - q_A)$$

# Cell 8
def combined_multiplicity(N_A, N_B, q_total):
    """Compute combined multiplicity for each energy split."""
    q_A_values = np.arange(0, q_total + 1)
    omega_combined = np.array([
        multiplicity(N_A, q_A) * multiplicity(N_B, q_total - q_A)
        for q_A in q_A_values
    ])
    return q_A_values, omega_combined

# Cell 9
# Small system: N_A = N_B = 10, q_total = 20
q_A, omega = combined_multiplicity(10, 10, 20)
probability = omega / omega.sum()

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].plot(q_A, omega)
axes[0].set_xlabel("Energy in solid A (q_A)")
axes[0].set_ylabel("Combined multiplicity")
axes[0].set_title("N_A = N_B = 10, q_total = 20")
axes[0].grid(True, alpha=0.3)

axes[1].plot(q_A, probability)
axes[1].set_xlabel("Energy in solid A (q_A)")
axes[1].set_ylabel("Probability P(q_A)")
axes[1].set_title("Probability distribution")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Most probable q_A: {q_A[np.argmax(omega)]}")
print(f"Peak probability: {probability.max():.4f}")

# Cell 10
# Large system: observe the sharpening
fig, ax = plt.subplots(figsize=(8, 5))

for N in [10, 50, 100, 200]:
    q_total = 2 * N  # Equal average energy per oscillator
    q_A_vals, omega_vals = combined_multiplicity(N, N, q_total)
    prob = omega_vals / omega_vals.sum()
    # Normalize x-axis to fraction of total energy
    ax.plot(q_A_vals / q_total, prob, label=f"N = {N}")

ax.set_xlabel("Fraction of energy in solid A (q_A / q_total)")
ax.set_ylabel("Probability")
ax.set_title("Sharpening of the multiplicity peak with system size")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Cell 11
# Quantify: standard deviation vs system size
sizes = [10, 20, 50, 100, 200, 500]
std_devs = []
relative_flucts = []

for N in sizes:
    q_total = 2 * N
    q_A_vals, omega_vals = combined_multiplicity(N, N, q_total)
    prob = omega_vals / omega_vals.sum()
    mean = np.sum(q_A_vals * prob)
    variance = np.sum((q_A_vals - mean)**2 * prob)
    std = np.sqrt(variance)
    std_devs.append(std)
    relative_flucts.append(std / mean)

plt.figure(figsize=(8, 5))
plt.loglog(sizes, relative_flucts, 'o-')
plt.loglog(sizes, 1.0 / np.sqrt(sizes), '--', label="1/sqrt(N) scaling")
plt.xlabel("System size N")
plt.ylabel("Relative fluctuation (sigma / mean)")
plt.title("Fluctuations vanish for large systems")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Cell 12 (Markdown)
# ## Step 4: Entropy and Temperature from Counting
#
# Entropy: $S = k \ln \Omega$ (using $k = 1$ natural units)
#
# Temperature: $\frac{1}{T} = \frac{\partial S}{\partial U}$ where $U = q \epsilon$

# Cell 13
def log_multiplicity(N, q):
    """Compute ln(Omega) using gammaln to avoid overflow for large systems."""
    return gammaln(q + N) - gammaln(q + 1) - gammaln(N)

# Verify against exact for small systems
print(f"Exact: ln(Omega(4,10)) = {log(multiplicity(4, 10)):.6f}")
print(f"Approx: ln(Omega(4,10)) = {log_multiplicity(4, 10):.6f}")

# Cell 14
# Entropy of two solids as a function of q_A
N_A, N_B, q_total = 100, 100, 200

q_A_vals = np.arange(0, q_total + 1)
S_A = np.array([log_multiplicity(N_A, q) for q in q_A_vals])
S_B = np.array([log_multiplicity(N_B, q_total - q) for q in q_A_vals])
S_total = S_A + S_B

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(q_A_vals, S_A, label="S_A (solid A)")
ax.plot(q_A_vals, S_B, label="S_B (solid B)")
ax.plot(q_A_vals, S_total, label="S_total = S_A + S_B", linewidth=2)
ax.axvline(q_A_vals[np.argmax(S_total)], color='red', linestyle='--',
           label=f"Equilibrium: q_A = {q_A_vals[np.argmax(S_total)]}")
ax.set_xlabel("Energy in solid A (q_A)")
ax.set_ylabel("Entropy S/k")
ax.set_title("Entropy maximization determines equilibrium")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Cell 15
# Derive temperature: 1/T = dS/dq (with epsilon = 1)
dS_A = np.diff(S_A)  # Numerical derivative
dS_B = -np.diff(S_B)  # Note: as q_A increases, q_B decreases

q_midpoints = q_A_vals[:-1] + 0.5

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(q_midpoints, dS_A, label="dS_A/dq_A = 1/T_A")
ax.plot(q_midpoints, dS_B, label="-dS_B/dq_A = 1/T_B")
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_xlabel("q_A")
ax.set_ylabel("1/T (inverse temperature)")
ax.set_title("Equilibrium: where temperatures are equal (curves cross)")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Find crossing point
crossing_idx = np.argmin(np.abs(dS_A - dS_B))
print(f"Temperatures equal at q_A ~ {q_midpoints[crossing_idx]:.1f}")
print(f"Maximum of S_total at q_A = {q_A_vals[np.argmax(S_total)]}")

You have now derived temperature from pure counting. Next up, you will use the Einstein model to predict a measurable quantity: heat capacity.

Plot the Einstein Heat Capacity

In the previous step, you proved that entropy maximization and equal temperatures are the same statement. Temperature emerged as a derivative of entropy with respect to energy.

Now you'll close the loop by computing a measurable physical quantity: heat capacity. This tells you how much energy a solid absorbs per degree of temperature increase, and the Einstein model makes a specific prediction that interpolates between quantum and classical behavior.

In this step, get ready to:
  • Compute the average energy per oscillator as a function of temperature.
  • Derive the Einstein heat capacity and identify its quantum and classical limits.
  • Compare the analytic formula to a numerical derivative to verify consistency.
Compute average energy as a function of temperature

So far you've worked with discrete energy quanta q and counted microstates directly. Now you'll connect that to continuous temperature.

For an Einstein solid in thermal equilibrium at temperature T, statistical mechanics gives a clean result for the average number of quanta per oscillator.

  • Add a new markdown cell below Cell 15 (your temperature-crossing output) to introduce the heat capacity section. Type the following into the markdown cell:
## Step 5: Heat Capacity of the Einstein Solid

For the Einstein solid at temperature $T$, the average number of quanta per oscillator is:

$$\langle q \rangle = \frac{1}{e^{\epsilon / kT} - 1}$$

The heat capacity per oscillator (in units of $k$) is:

$$\frac{C_V}{Nk} = \left(\frac{\epsilon}{kT}\right)^2 \frac{e^{\epsilon/kT}}{(e^{\epsilon/kT} - 1)^2}$$

What do these formulas describe?

The energy formula is identical to the Planck distribution for photons or phonons. Each oscillator holds an average of <q> quanta, which goes from 0 at low temperature (all oscillators frozen in the ground state) to roughly kT/epsilon at high temperature (classical behavior).

The heat capacity formula is the temperature derivative of the energy. It captures how responsive the solid is to temperature changes.

  • Add a new code cell below the markdown cell and define the physics functions by entering this code:
def einstein_heat_capacity(T, theta_E=1.0):
    """Dimensionless heat capacity C_V / (Nk) for the Einstein solid.
    
    T and theta_E in same units. theta_E = epsilon/k is the Einstein temperature.
    """
    x = theta_E / T
    return x**2 * np.exp(x) / (np.exp(x) - 1)**2

def einstein_energy(T, theta_E=1.0):
    """Average energy per oscillator in units of epsilon."""
    x = theta_E / T
    return 1.0 / (np.exp(x) - 1)

What does this code do?

The einstein_energy function computes the Bose-Einstein occupation number for a single oscillator at temperature T. The variable x = theta_E / T is the dimensionless ratio that controls the physics.

The einstein_heat_capacity function returns the analytic derivative dU/dT in dimensionless form. Both functions use theta_E = 1.0 as the default Einstein temperature, so T is measured in units of theta_E.

  • Run Cell 17. No output appears because you are only defining functions.

Seeing a NameError for np?

Make sure you ran Cell 2 (the imports cell) earlier in this session. If you restarted your kernel, re-run all cells from the top by clicking **Run**, then **Run All Above Selected Cell**.

help me fix a NameError in my Jupyter notebook

Understand the two temperature limits

The Einstein heat capacity makes a prediction that was historically groundbreaking. It explains why solids stop absorbing heat at low temperatures, something classical physics could not account for.

What are the two limits?

At high temperature (T much greater than theta_E), each oscillator has many quanta of energy and behaves classically. The heat capacity approaches 1 in units of Nk per oscillator. This is the Dulong-Petit limit, observed experimentally in most solids at room temperature.

At low temperature (T much less than theta_E), oscillators are frozen in their ground state because thermal energy is too small to excite even one quantum. The heat capacity drops to zero exponentially. This is quantum freezeout.

Compare numerical derivative to analytic formula

The ultimate check is to compute C_V two ways: from the analytic formula and from a numerical derivative of the energy curve. If the functions are correct, both methods produce identical curves.

  • Add a new code cell below Cell 17 and compute the energy, heat capacity, and generate the two-panel plot by entering this code:
# Plot heat capacity
T_values = np.linspace(0.01, 5.0, 500)  # In units of theta_E
C_analytic = einstein_heat_capacity(T_values, theta_E=1.0)

# Numerical derivative for comparison
U_values = einstein_energy(T_values, theta_E=1.0)
C_numerical = np.gradient(U_values, T_values)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Energy plot
axes[0].plot(T_values, U_values)
axes[0].set_xlabel("T / theta_E")
axes[0].set_ylabel("Average energy per oscillator (units of epsilon)")
axes[0].set_title("Energy vs Temperature")
axes[0].axhline(1, color='gray', linestyle='--', alpha=0.5, label="High-T: <q> ~ T")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Heat capacity plot
axes[1].plot(T_values, C_analytic, label="Analytic Einstein formula", linewidth=2)
axes[1].plot(T_values, C_numerical, '--', label="Numerical dU/dT", alpha=0.7)
axes[1].axhline(1, color='red', linestyle=':', label="Dulong-Petit limit (C=Nk)")
axes[1].set_xlabel("T / theta_E")
axes[1].set_ylabel("C_V / (Nk)")
axes[1].set_title("Einstein Heat Capacity")
axes[1].set_ylim(0, 1.1)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("Low T behavior: C_V -> 0 exponentially (quantum freezeout)")
print(f"C_V at T=0.1: {einstein_heat_capacity(0.1):.6f}")
print(f"C_V at T=5.0: {einstein_heat_capacity(5.0):.6f} (approaching 1.0)")

What does this code do?

The code computes heat capacity two independent ways and plots them on the same axes.

  • C_analytic uses the closed-form Einstein formula directly.
  • C_numerical uses np.gradient to take a finite-difference derivative of the energy curve with respect to temperature.
  • The red dotted line at C_V/(Nk) = 1 marks the classical Dulong-Petit limit that all solids approach at high temperature.

If the two curves overlap perfectly, the energy function and heat capacity function are consistent with each other.

  • Run Cell 18.
  • You should see a two-panel figure. The left panel shows energy rising from zero and approaching linearity at high T. The right panel shows both heat capacity curves sitting directly on top of each other, rising from near zero toward the Dulong-Petit limit at 1.
  • Check the printed output below the plot. C_V at T=0.1 should be extremely small (near zero), and C_V at T=5.0 should be close to 1.0.

Curves don't overlap?

If the dashed numerical curve diverges from the solid analytic curve, check that T_values starts at 0.01 (not zero, which causes division by zero). Also verify you passed T_values as the second argument to np.gradient so it knows the spacing between points.

help me debug my Einstein heat capacity plot

Why does this matter?

You just reproduced Einstein's 1907 result that resolved the classical heat capacity puzzle. Classical physics predicted that all solids should have the same heat capacity at all temperatures. Experiments showed heat capacity drops toward zero at low temperatures, and Einstein's quantum model was the first to explain why.

✔️ Awesome, I've got everything!

Great. Double check you have saved your notebook before moving on.

ⓧ I'd like to double check the full code

Here is the complete notebook up to and including this step. Compare your cells against this:

# Cell 1 (Markdown)
# # Entropy from Counting: Two Einstein Solids in Thermal Contact
# 
# Following Schroeder's *Thermal Physics* Chapter 2 approach:
# start by counting microstates, then watch thermodynamics emerge.

# Cell 2
import numpy as np
import matplotlib.pyplot as plt
from math import comb, log
from scipy.special import gammaln

# Cell 3 (Markdown)
# ## Step 2: Multiplicity of a Single Einstein Solid
#
# An Einstein solid has N quantum harmonic oscillators sharing q indistinguishable
# energy quanta. The multiplicity (number of microstates) is the "stars and bars"
# combinatorial formula:
#
# $$\Omega(N, q) = \binom{q + N - 1}{q}$$

# Cell 4
def multiplicity(N, q):
    """Number of microstates for N oscillators sharing q energy quanta."""
    return comb(q + N - 1, q)

# Verify against Schroeder's examples
print(f"Omega(4, 10) = {multiplicity(4, 10)}")  # Should be 286
print(f"Omega(4, 2) = {multiplicity(4, 2)}")    # Should be 10
print(f"Omega(3, 3) = {multiplicity(3, 3)}")    # Should be 10

# Cell 5
# Explore small systems
for N in [3, 4, 10]:
    for q in [2, 5, 10]:
        print(f"N={N}, q={q}: Omega = {multiplicity(N, q)}")

# Cell 6
# Plot multiplicity growth for fixed N
N = 50
q_values = np.arange(0, 101)
omega_values = [multiplicity(N, q) for q in q_values]

plt.figure(figsize=(8, 5))
plt.semilogy(q_values, omega_values)
plt.xlabel("Energy quanta q")
plt.ylabel("Multiplicity (log scale)")
plt.title(f"Multiplicity of an Einstein Solid (N={N} oscillators)")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Cell 7 (Markdown)
# ## Step 3: Two Solids in Thermal Contact
#
# Two solids A and B share a fixed total energy $q_{total}$.
# The combined multiplicity for a given split is:
#
# $$\Omega_{total}(q_A) = \Omega_A(q_A) \times \Omega_B(q_{total} - q_A)$$

# Cell 8
def combined_multiplicity(N_A, N_B, q_total):
    """Compute combined multiplicity for each energy split."""
    q_A_values = np.arange(0, q_total + 1)
    omega_combined = np.array([
        multiplicity(N_A, q_A) * multiplicity(N_B, q_total - q_A)
        for q_A in q_A_values
    ])
    return q_A_values, omega_combined

# Cell 9
# Small system: N_A = N_B = 10, q_total = 20
q_A, omega = combined_multiplicity(10, 10, 20)
probability = omega / omega.sum()

fig, axes = plt.subplots(1, 2, figsize=(12, 4))

axes[0].plot(q_A, omega)
axes[0].set_xlabel("Energy in solid A (q_A)")
axes[0].set_ylabel("Combined multiplicity")
axes[0].set_title("N_A = N_B = 10, q_total = 20")
axes[0].grid(True, alpha=0.3)

axes[1].plot(q_A, probability)
axes[1].set_xlabel("Energy in solid A (q_A)")
axes[1].set_ylabel("Probability P(q_A)")
axes[1].set_title("Probability distribution")
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print(f"Most probable q_A: {q_A[np.argmax(omega)]}")
print(f"Peak probability: {probability.max():.4f}")

# Cell 10
# Large system: observe the sharpening
fig, ax = plt.subplots(figsize=(8, 5))

for N in [10, 50, 100, 200]:
    q_total = 2 * N  # Equal average energy per oscillator
    q_A_vals, omega_vals = combined_multiplicity(N, N, q_total)
    prob = omega_vals / omega_vals.sum()
    # Normalize x-axis to fraction of total energy
    ax.plot(q_A_vals / q_total, prob, label=f"N = {N}")

ax.set_xlabel("Fraction of energy in solid A (q_A / q_total)")
ax.set_ylabel("Probability")
ax.set_title("Sharpening of the multiplicity peak with system size")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Cell 11
# Quantify: standard deviation vs system size
sizes = [10, 20, 50, 100, 200, 500]
std_devs = []
relative_flucts = []

for N in sizes:
    q_total = 2 * N
    q_A_vals, omega_vals = combined_multiplicity(N, N, q_total)
    prob = omega_vals / omega_vals.sum()
    mean = np.sum(q_A_vals * prob)
    variance = np.sum((q_A_vals - mean)**2 * prob)
    std = np.sqrt(variance)
    std_devs.append(std)
    relative_flucts.append(std / mean)

plt.figure(figsize=(8, 5))
plt.loglog(sizes, relative_flucts, 'o-')
plt.loglog(sizes, 1.0 / np.sqrt(sizes), '--', label="1/sqrt(N) scaling")
plt.xlabel("System size N")
plt.ylabel("Relative fluctuation (sigma / mean)")
plt.title("Fluctuations vanish for large systems")
plt.legend()
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Cell 12 (Markdown)
# ## Step 4: Entropy and Temperature from Counting
#
# Entropy: $S = k \ln \Omega$ (using $k = 1$ natural units)
#
# Temperature: $\frac{1}{T} = \frac{\partial S}{\partial U}$ where $U = q \epsilon$

# Cell 13
def log_multiplicity(N, q):
    """Compute ln(Omega) using gammaln to avoid overflow for large systems."""
    return gammaln(q + N) - gammaln(q + 1) - gammaln(N)

# Verify against exact for small systems
print(f"Exact: ln(Omega(4,10)) = {log(multiplicity(4, 10)):.6f}")
print(f"Approx: ln(Omega(4,10)) = {log_multiplicity(4, 10):.6f}")

# Cell 14
# Entropy of two solids as a function of q_A
N_A, N_B, q_total = 100, 100, 200

q_A_vals = np.arange(0, q_total + 1)
S_A = np.array([log_multiplicity(N_A, q) for q in q_A_vals])
S_B = np.array([log_multiplicity(N_B, q_total - q) for q in q_A_vals])
S_total = S_A + S_B

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(q_A_vals, S_A, label="S_A (solid A)")
ax.plot(q_A_vals, S_B, label="S_B (solid B)")
ax.plot(q_A_vals, S_total, label="S_total = S_A + S_B", linewidth=2)
ax.axvline(q_A_vals[np.argmax(S_total)], color='red', linestyle='--',
           label=f"Equilibrium: q_A = {q_A_vals[np.argmax(S_total)]}")
ax.set_xlabel("Energy in solid A (q_A)")
ax.set_ylabel("Entropy S/k")
ax.set_title("Entropy maximization determines equilibrium")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Cell 15
# Derive temperature: 1/T = dS/dq (with epsilon = 1)
dS_A = np.diff(S_A)  # Numerical derivative
dS_B = -np.diff(S_B)  # Note: as q_A increases, q_B decreases

q_midpoints = q_A_vals[:-1] + 0.5

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(q_midpoints, dS_A, label="dS_A/dq_A = 1/T_A")
ax.plot(q_midpoints, dS_B, label="-dS_B/dq_A = 1/T_B")
ax.axhline(0, color='gray', linestyle='-', alpha=0.3)
ax.set_xlabel("q_A")
ax.set_ylabel("1/T (inverse temperature)")
ax.set_title("Equilibrium: where temperatures are equal (curves cross)")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# Find crossing point
crossing_idx = np.argmin(np.abs(dS_A - dS_B))
print(f"Temperatures equal at q_A ~ {q_midpoints[crossing_idx]:.1f}")
print(f"Maximum of S_total at q_A = {q_A_vals[np.argmax(S_total)]}")

# Cell 16 (Markdown)
# ## Step 5: Heat Capacity of the Einstein Solid
#
# For the Einstein solid at temperature $T$, the average number of quanta per oscillator is:
#
# $$\langle q \rangle = \frac{1}{e^{\epsilon / kT} - 1}$$
#
# The heat capacity per oscillator (in units of $k$) is:
#
# $$\frac{C_V}{Nk} = \left(\frac{\epsilon}{kT}\right)^2 \frac{e^{\epsilon/kT}}{(e^{\epsilon/kT} - 1)^2}$$

# Cell 17
def einstein_heat_capacity(T, theta_E=1.0):
    """Dimensionless heat capacity C_V / (Nk) for the Einstein solid.
    
    T and theta_E in same units. theta_E = epsilon/k is the Einstein temperature.
    """
    x = theta_E / T
    return x**2 * np.exp(x) / (np.exp(x) - 1)**2

def einstein_energy(T, theta_E=1.0):
    """Average energy per oscillator in units of epsilon."""
    x = theta_E / T
    return 1.0 / (np.exp(x) - 1)

# Cell 18
# Plot heat capacity
T_values = np.linspace(0.01, 5.0, 500)  # In units of theta_E
C_analytic = einstein_heat_capacity(T_values, theta_E=1.0)

# Numerical derivative for comparison
U_values = einstein_energy(T_values, theta_E=1.0)
C_numerical = np.gradient(U_values, T_values)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Energy plot
axes[0].plot(T_values, U_values)
axes[0].set_xlabel("T / theta_E")
axes[0].set_ylabel("Average energy per oscillator (units of epsilon)")
axes[0].set_title("Energy vs Temperature")
axes[0].axhline(1, color='gray', linestyle='--', alpha=0.5, label="High-T: <q> ~ T")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# Heat capacity plot
axes[1].plot(T_values, C_analytic, label="Analytic Einstein formula", linewidth=2)
axes[1].plot(T_values, C_numerical, '--', label="Numerical dU/dT", alpha=0.7)
axes[1].axhline(1, color='red', linestyle=':', label="Dulong-Petit limit (C=Nk)")
axes[1].set_xlabel("T / theta_E")
axes[1].set_ylabel("C_V / (Nk)")
axes[1].set_title("Einstein Heat Capacity")
axes[1].set_ylim(0, 1.1)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

print("Low T behavior: C_V -> 0 exponentially (quantum freezeout)")
print(f"C_V at T=0.1: {einstein_heat_capacity(0.1):.6f}")
print(f"C_V at T=5.0: {einstein_heat_capacity(5.0):.6f} (approaching 1.0)")

Secret mission

The Boltzmann Distribution from Random Exchanges

Everything so far has been exact enumeration. What if you just let oscillators swap energy at random? In this mission, you will simulate a single Einstein solid reaching thermal equilibrium through pure randomness and watch the Boltzmann distribution emerge from nothing but statistics.

Clean Up Your Resources

Clean Up Your Resources

This project runs entirely locally, so there are no ongoing costs. Decide whether to keep your resources running, pause them to come back later, or delete them entirely.

Resources you used:

  • JupyterLab server (currently running in your terminal).
  • einstein_solids.ipynb notebook file.
  • Python packages installed via pip: jupyterlab, numpy, matplotlib, scipy.

Keep everything running

No action needed. Choose this if you're still actively building or want to revisit the notebook later today.

  • Your JupyterLab server stays running in the terminal where you launched it.
  • Your notebook and all installed packages remain available.
  • To reopen the notebook later, navigate to the JupyterLab tab in your browser or run jupyter lab again from the same directory.

Pause - I'll come back to this later

Shut down the JupyterLab server to free up memory. Your notebook file and packages stay installed for next time.

  • Switch to the terminal where JupyterLab is running.
  • Press Ctrl+C twice to stop the server.
  • Your einstein_solids.ipynb file and all pip packages remain on disk.
  • To resume later, open your terminal, navigate to the folder containing the notebook, and run jupyter lab.

Delete - I don't want to use this again

Remove all project resources and uninstall the packages you added.

  • Switch to the terminal where JupyterLab is running and press Ctrl+C twice to stop the server.
  • Delete the notebook file by running:
rm einstein_solids.ipynb
  • Uninstall the Python packages you installed for this project by running:
pip uninstall jupyterlab numpy matplotlib scipy -y

What about Python itself?

Python is a general-purpose tool you will likely use for other projects. These instructions only remove the packages installed specifically for this project, not Python itself.

Nice Work!

Nice Work!

Nice work! You just built a computational proof that the second law of thermodynamics is not a separate axiom but a statistical certainty.

You've learned how to:

  • Demonstrate that the second law of thermodynamics emerges from pure combinatorics by computing multiplicity for two interacting Einstein solids and showing the peak becomes overwhelmingly sharp for large systems.
  • Derive entropy and temperature from microstate counting, proving that thermal equilibrium is simply entropy maximization.
  • Compute the Einstein heat capacity curve showing both quantum freezeout at low temperatures and the classical Dulong-Petit limit at high temperatures.
  • Secret Mission: Simulated random energy exchanges between oscillators and watched the Boltzmann distribution emerge from nothing but statistics.

Ready to quiz yourself?