The Linear Molecule
We begin our spectroscopic adventure with the most simple case, a linear molecule. Here, we will use OCS as an example as it has a considerable dipole moment and is often used as a guinea pig molecule in the laboratory.
Theory
The most simple assumption is that of a rigid linear molecule. As a linear molecule has no rotation around the $a$-axis and $I_b = I_c = I$ its rotational Hamiltonian is simplified to $$ \mathbf{H}_\text{rot} = \frac{\mathbf{J}^2}{I} $$ and the rotational energies are its eigenvalues $$ E_\text{rot}/h = BJ(J+1) $$ as shown on the left hand side of the figure. For one-photon transitions, the change in angular momentum is limited to the photon’s angular momentum of $\pm \hbar$ resulting in $\Delta J = \pm 1$. The resulting spectrum consists of equidistant transitions $$ \nu_{J+1 \leftarrow J} = 2B(J+1) $$ with a distance of $2B$ as shown on the right hand side of the figure.
When adding centrifugal distortion for a non-rigid molecule, the energy expression is expanded into $$ E_\text{rot}/h =BJ(J+1)-D\left(J(J+1)\right)^{2} + \text{...} $$ with the centrifugal distortion constant $D$. For linear molecules, $D$ is positive as the internuclear distance is increased by the centrifugal forces. As a result, the effective rotational constant $B_{eff}=B-DJ(J+1)$ and the transition frequencies $$ \nu_{J+1 \leftarrow J} = 2B(J+1) - 4 D (J+1)^3 + \text{...} $$ are decreased in comparison to the rigid rotor.
The main results are, that rotational spectra of linear molecules have almost equidistant patterns (as $D$ is a few magnitudes smaller than $B$) and that neglecting higher centrifugal distortion constants results in transition frequencies being predicted to slightly higher frequencies.
Approximation of the Rotational Constant $B$
Before we look at a spectrum, we will approximate the rotational constant $B$ of OCS. As OCS is a linear molecule, its geometry is defined by two distances, the CO distance $r_\text{CO}$ and the CS distance $r_\text{CS}$. They can be approximated via the rotational constants of the two diatomic molecules CO and CS ($B_\text{CO} = 57635.968 \text{ MHz}$ and $B_\text{CS} = 24495.560 \text{ MHz}$). Please calculate $r_\text{CO}$ from $B_\text{CO}$ and $r_\text{CS}$ and $B_\text{CS}$ from $r_\text{CS}$. Then calculate the rotational constant of OCS, assuming that $r_\text{CO}$ and $r_{CS}$ would not change when going from the two diatomic molecules to OCS (which of course is not true but allows us for a very simple estimation).
Solution
For a two atomic molecule with one atom with mass $m_1$ at the origin and the other atom with mass $m_2$ at distance $r$, the center of mass $c$ is $$ c = \frac{r \cdot m_2}{m_1 + m_2} $$ Then the moment of inertia is given by $$ I = \sum_i m_i r_i^2 = \mu \cdot r^2 $$ With the reduced mass $\mu = \frac{m_1 \cdot m_2}{m_1+m_2}$.
Furthermore, the rotational constant in units of a frequency is given as $$ B = \frac{h}{8 \pi^2 I} $$ Thus, the radius $r$ for such a two atomic molecule is given as $$ r = \sqrt{\frac{h}{8 \pi^2 \mu B}} $$
With the masses of $m_\text{C} = 12$, $m_\text{0} = 16$, and $m_\text{S} = 32$ the bond distances for CO and CS can be calculated as $r_\text{CO} = 1.1308 \cdot 10^{-10} \text{ m}$ and $r_\text{CS} = 1.5375 \cdot 10^{-10} \text{ m}$. From there we can estimate the rotational constant of OCS via the moment of inertia as $B=6318.8379 \text{ MHz}$.
The Python implementation of the solution is given below.
Solution
import numpy as np
B_CO = 57635.968E6
B_CS = 24495.560E6
h = 6.626e-34
m0 = 1.66054e-27
ms = {'C': 12, 'O': 16, 'S': 32}
mu_CO = (ms['C'] * ms['O']) / (ms['C'] + ms['O']) * m0
mu_CS = (ms['C'] * ms['S']) / (ms['C'] + ms['S']) * m0
r_CO = np.sqrt(h/(8 * np.pi**2 * mu_CO * B_CO))
r_CS = np.sqrt(h/(8 * np.pi**2 * mu_CS * B_CS))
print(f'The distance for CO is {r_CO} m')
print(f'The distance for CS is {r_CS} m')
masses = np.array((ms['O'], ms['C'], ms['S']))
positions = np.array((0, r_CO, r_CO + r_CS))
c = np.sum(masses * positions) / np.sum(masses)
positions -= c
I = np.sum(masses * m0 * positions**2)
B = h / (8 * np.pi**2 * I)
print(f'The rotational constant for OCS is about {B / 1E6:.4f} MHz')
Exploring the spectrum
For this exercise, we will use a synthetically generated spectrum. This has a few didactic advantages. Download the spectrum here. It is in *.csv format with the frequencies in MHz in the first column and the intensities in arbitrary units in the second column.
Open LLWP by opening a command line window and running
llwp
Then load the spectrum into LLWP either via drag-and-drop or via the Menu (Files > Add Files).
Navigating the Spectrum
Then set the number of plots to one (Ctrl + N then enter 1 in the dialog), the center of the plot to 100 GHz (Ctrl + G then enter 100000 in the dialog), and the width of the plot to 200 GHz (Ctrl + W then enter 2E5 in the dialog). The plot should look now as shown in the figure.
Now you can navigate the spectrum by various means: - use the in, out, left, right buttons above the plot (or W, S, A, D on your keyboard, respectively) - set different width and center positions as described above - Hold shift and select a region with your mouse to zoom into that region - Use the scroll wheel of your mouse (while hovering over the plot) to zoom in and out
Loomis-Wood Plots
To see how well our estimated $B$ value compares to the experimental spectrum, we will now go to a Loomis-Wood representation of the spectrum. The idea is that same-width excerpts of the spectrum around the predicted positions of a transition series are aligned vertically.
Here, we will look at the five transitions $J = 1 \leftarrow 0$ to $J=5 \leftarrow 4$. To do so, perform the following steps: - Change the number of plots to 5 (Ctrl + N then enter 5 in the dialog) - Change the relative offset of the plots to 0 (Ctrl + G then enter 0 in the dialog) - Change the width of each plot to 10 GHz (Ctrl + W then enter 10E3 in the dialog) - Go to the Reference Series window and select on top Expression - Enter the following expression: 2 * 6318.8379 * (N+1+N0) - Press Update - Change the annotations of the plots - Open the config (Shift + 0 or View > Config), increase the size of the window if required - In the search field, search for plot_annotationfstring and change its value to {x:.2f} - This will display the reference position of each plot in the top right corner
You should see now a plot similar to the following figure.
You see five excerpts of the spectrum. They are centered around the predicted positions of the $J = 1 \leftarrow 0$ to $J=5 \leftarrow 4$ transitions (from bottom to top). The reference positions (the centers of the five plots) are determined by the expression we entered into the Reference Series window.
The experimental spectrum deviates toward lower frequencies (to the left) when going from low to high $J$ values (from bottom to top). This clearly indicates, that our estimated $B$ value is too high.
Fitting the Center Position of a Line and Improving the Model
To estimate a better $B$ value, we will determine the position of the lowest transition and divide its value by 2. - Zoom into the strong line of the lowest plot to see its profile (either with in, out, left, right or by holding shift and selecting an area around the line) - Choose the right fit function, here a simple Gaussian (Fit > Choose Fit Function > Gauss) - Select the area around the line with the mouse (without holding shift)
You should see a few things appear: - the fitted lineshape on top of the experimental spectrum - the center position on top of the experimental spectrum indicated by a dashed line - a new entry in the New Assignments window
Press Resize in the New Assignments window to resize the columns of the table. You can see the center position in the Freq column but it will also be copied to your clipboard by default. Divide the value by two
B_exp = 12162.982151012036 / 2
print(B_exp)
Now update the expression in the Reference Series window accordingly
2 * 6081.49 * (N+1+N0)
and press Update. Set the relative offset of the plots to zero (Ctrl + G then enter 0 in the dialog) and their back to 10 GHz (Ctrl + W then enter 10E3 in the dialog). As you can see, the new $B$ value describes the series much better and the series almost form a vertical trend. Lets zoom into this Loomis-Wood plot to see that the description still is not perfect (a value of 5 MHz for the width is appropriate). The maximum deviations are now a lot smaller (about 1 MHz) and the trend is less of a straight line and more of a curve. To see this more pronounce, go to 10 plots and a width of 15 MHz.
The deviations are a result of centrifugal distortion which we have omitted in our very simple model. As described in the theoretical section, centrifugal distortion is included as higher order terms. You can try to estimate $D$ using the equation $$ \nu_{J+1 \leftarrow J} = 2B(J+1) - 4 D (J+1)^3 + \text{...} $$ and the center frequencies for the transitions $J = 1 \leftarrow 0$ and $J=10 \leftarrow 9$. This should yield $D=1.29 \text{ kHz}$. Update the expression in the Reference Series window accordingly
2 * 6081.49 * (N+1+N0) - 4 * 1.29E-03 * (N+1+N0)**3
and press update to see how well this model reproduces the first ten transitions.
Solution
v_exp_0 = 12162.982151012036
v_exp_9 = 121624.63819944281
B_exp = 12162.982151012036 / 2
v_calc_9 = 2 * B_exp * 10
D_exp = (v_exp_9 - v_calc_9) / (4 * (9 + 1)**3)
print(f'{B_exp=:.2f}')
print(f'{D_exp=:.2e}')
Modelling OCS in SPFIT
We have seen that we can calculate the rotational parameter $B$ and the centrifugal distortion parameter $D$ quite easily and quite accurately from the experimental spectrum. However, ideally we would assign all visible transitions to the correct quantum numbers and then fit the Hamiltonian to these assignments and to create new predictions based on this updated Hamiltonian. We will learn how to perform these tasks with SPFIT and SPCAT for a simple linear molecule in the next section.