In the outer layers of giant planets, molecular hydrogen occurs in two different states: para- and ortho-hydrogen: https://en.wikipedia.org/wiki/Spin_isomers_of_hydrogen Both have different energies, which means their fractions are temperature dependent. Spacecrafts can measure their fractions remotely. If a deviation from thermodynamic equilibrium is detected, it implies that convective forces brought up hydrogen from deeper, hotter layers at a faster rate than it could equilibrate. Therefore such measurements are conducted to gauge the speed of convection in giant planet atmospheres.
In this assignment, we only want to compute the fraction of para-hydrogen as function of temperature in thermodynamic equilibrium. Please review the lecture of B. Militzer for further information.
We assume all electrons are in the ground state because the energies of the excited state are too high for the temperature range under consideration. According to the National Instiute of Standards (NIST) http://webbook.nist.gov/cgi/cbook.cgi?ID=C1333740&Units=SI&Mask=1000#Diatomic the energy of eigenstates with vibrational quantum number $m$ and rotational quantum number $J$ is given by $$ \epsilon_{m,J} = \omega_e (m+\tfrac{1}{2}) - \omega_e x_e \left(m+\tfrac{1}{2}\right)^2 + B_e J(J+1) - D_e J^2(J+1)^2 - \alpha_e (m+\tfrac{1}{2}) J(J+1)\,.$$ In typical spectroscopy notation, the constants are given in wave numbers with unit of cm$^{-1}$: $W_e = 4401.213, W_eX_e = 121.336, B_e = 60.853, D_e = 0.0471$ and $ \alpha_e = 3.0622\,.$
The corresponding partition function in the canonical ensemble is given by $$ Z(T) = \sum_{m,J} g_{NS} \times g_J \times e^{-\epsilon_{m,J}/k_BT},$$ which is a summation of the statistical weights of all eigenstates with $m \ge 0$ and $J \ge 0$. The fraction of hydrogen molecules in the para state is equal to the fraction of states with even $J$ because of these properties:
Para hydrogen: ground state, antisymmetric nuclear spin wf, symmetric rotational wf, even $J$ only
Ortho hydrogen: excited state, symmetric nuclear spin wf, antisymmetric rotational wf, odd $J$ only
$g_{NS}$ represents the degeneracy of the nuclear spin state. For hydrogen, it is 3 for odd $J$ and 1 for even $J$.
$g_J = 2J+1$ represents the degeneracy of the rotational state.
Although it is not needed for this assignment, it might be useful to know that the free energy, $F$, can be derived directly from the partition function $$e^{-\beta F} = Z$$ and that the internal energy, $E$, is given by the following statistical average of all eigen energies, $$ E = \left<\epsilon_{m,J}\right> = \frac{1}{Z} \sum_{m,J} \epsilon_{m,J} \times g_{NS} \times g_J \times e^{-\epsilon_{m,J}/k_BT},$$
import numpy as np
import matplotlib.pyplot as plt
eConst = 1.602176462e-19
mu0 = 4.0*np.pi*1e-7
c = 299792458.0 # m/s
eps0 = 1.0/(mu0*c*c)
fConst = eConst*eConst/(4.0*np.pi*eps0)
h = 6.62606876e-34
hBar = h/(2.0*np.pi)
kb = 1.3806503e-23
u = 1.66053873e-27 # atomic mass unit in kg
me = 9.10938188e-31
mp = 1.67262158e-27
md = 1.99900750083*mp
a0 = hBar*hBar/(fConst*me)
Ha = fConst*fConst*me/hBar/hBar
Ry = Ha/2.0
NA = 6.0221367e23
A = 1e-10 # Angstroem
StephanBoltzmannConstant = 5.670367e-8
# For a given temperature in Kelvin, this function should return
# the fraction of hydrogen molecules that are in a para state
# Replace all ... by meaning Python commands. Introduce new
# variables as needed.
def ParaFraction(T):
#print('T=',T,'beta=',beta)
z = 0.0
zPara = 0.0
u = 0.0
zLast = 0.0
mMax = 18
# loop over vibrational states. Remember to break out early of this loop below
for m in range(mMax):
epsilonLast = 0.0
jMax = 100
# loop over vibrational states. Remember to break out early of this loop below
for j in range(jMax):
epsilon = # For given m and j, compute energy of eigen state
#print(m,j,epsilon)
# ironically, the equation for epsilon(m,j) becomes inaccurate for large j
# So, if you detect that epsilon became less the value for the previous j, break out of
# this loop without considering the new epsilon(m,j) value.
if (...):
break
# For given m and j, compute this state contribution to the partition function
ys = ...
# add it to variable z
z += ys
# if this state is a para state, add it to zPara
if (...):
zPara += ys
# this marks the end of the 'j' loop. It is time to relax. Take a deep breath.
# After completing the 'j' loop, 'z' has a slightly larger value than before.
# If m>5 and the changes in the 'z' value was less than 1e-10, break of out of the 'm' loop now
if (m>5 and ...): break;
paraFraction = zPara / z
print('-->',T,paraFraction)
return paraFraction
Assume your function works correctly, please answere these four questions:
(1) What is the fraction of para hydrogen in the low-temperature limit:
...
(2) What is the fraction of para hydrogen in the high-temperature limit:
...
(3) For which temperature does the fraction of para hydrogen assume a value of 1/2? (accuracy 5 K ok)
...
(4) Plot the para fraction as function of temperature. Choose a temperature range so that one can see the system approach the low and high temperature limits.