2.2. Boltzmann Distribution Harmonic Oscillator

In this part we will create a simple computer program to compute the Boltzmann distribution of a fictious harmonic oscillator.

Exercise 5

Modify the code below to calculate the occupancy of each state within the harmonic oscillator system. Present the entire code file within your report and comment upon the main features using #. Consider different (reduced) temperatures, 0.5, 1, 2 and 3, and energy levels set up to 10 (recall \(\epsilon =\) integer values ). You will implement different StateOccupancy() functions, correspondng to different degeneracies. Note that in the degeneracy cases like \(s+1\), \(s\) is the index of the energy level, and its degeneracy is \((s + X)\), where \(X\) assumes given values in this Exercise.

Exercise 6

What do you note as the temperature is increased in terms of state occupancy? Does it make sense physically?

Exercise 7

Compare different degeneracies (no degeneracy, \(s+1\), and \(s+2\)) for different reduced temperatures. What can you infer?

import numpy as np
import matplotlib.pyplot as plt
import sys
sys.path.append("..")
import helpers
helpers.set_style()

plt.style.use(['seaborn-ticks', 'seaborn-talk'])

# MODIFY HERE
# set number of energy levels and temperatures here
n_energy_levels= 0
reducedTemperatures=[0, 0, 0, 0]


def calculateStateOccupancy(T, i):
    # No degeneracy
    return 1  # MODIFY HERE

def calculateStateOccupancy_s1(T, i):
    # Degeneracy s+1
    return 1  # MODIFY HERE

def calculateStateOccupancy_s2(T, i):
    # Degeneracy s+2
    return 1  # MODIFY HERE


functions ={
    "no_degeneracy": calculateStateOccupancy,
    "s+1": calculateStateOccupancy_s1,
    "s+2": calculateStateOccupancy_s2
}

for f in functions.keys():
    fig, ax =plt.subplots(1)
    ax.set_title(f)
    ax.set_xlabel("Energy level")
    ax.set_ylabel("Occupancy")

    calculateOccupancy=functions[f]

    for reducedTemperature in reducedTemperatures:
        distribution = [] # For each state there is one entry
        partition_function=np.float64(0.0)

        for i in range(n_energy_levels):
            stateOccupancy=calculateOccupancy(reducedTemperature, i)
            distribution.append(1.0) # MODIFY HERE
            partition_function+=1.0  # MODIFY HERE

        ax.plot(distribution/partition_function, label=reducedTemperature)

    ax.legend()
    plt.show()

Exercise 8

Modify your program such that the state occupancy and partition function are calculated for a linear rotor with moment of inertia \(I\) (copy your code in the cell below, and adapt it for the linear rotor case. Present also this code within your report). Compare your results at different temperatures with the approximate result (2.25).

Hint: considering more than four reduced temperatures can be a good idea, you can easily obtain evenly spaced numbers over a specified interval using the numpy function np.linspace(start, stop, number_of_points). Note that if the higher temperatures considered are far from the highest temperature which was considered so far, more energy levels may become accessible to the system, and you should keep them into account.

Is this an approximated result for which temperature limit?

(2.25)\[ Z = \frac{2I}{\beta\hbar^2} \]

using

\[ \frac{I}{\hbar^2} = 1. \]

Note that the energy levels of a linear rotor are:

\[ U = J(J+1)\frac{\hbar^2}{2I} \]

with \(J= 0, 1, 2, \dots, \infty\) and the degeneracy of level J is \(2J +1\).

# INSERT YOUR WORKING CODE HERE
# ADAPT it for the linear rotor case 
# defining a new calculateStateOccupancy_rotor function and a "linear rotor" key

# Then run the following cell (after having edited it where indicated!)

for f in functions.keys():
    fig, ax =plt.subplots(1)
    ax.set_title(f)
    ax.set_xlabel("Reduced Temperature")
    ax.set_ylabel("Partition Function")

    calculateOccupancy=functions[f]
    Z = []  # This variable will store the partition function 
            # for different reduced temperatures
    
    for reducedTemperature in reducedTemperatures:
        distribution = [] # For each state there is one entry
        partition_function=np.float64(0.0)

        for i in range(n_energy_levels):
            stateOccupancy=calculateOccupancy(reducedTemperature, i)
            distribution.append(1.0) # MODIFY HERE
            partition_function+=1.0  # MODIFY HERE
        Z.append(1.0) # MODIFY HERE
    ax.plot(reducedTemperatures, Z, label="Z Exact")
    ax.plot(reducedTemperatures, Z, label="Z Approximated") # MODIFY HERE
                                                            # insert approximated Z
    
    ax.legend()
    plt.show()