The Linear Molecule

We will begin our spectroscopic adventure with the simplest case: a linear molecule. For this example, we will use OCS, which has a significant dipole moment and is frequently employed as a test molecule in laboratories.

Theory

The simplest 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}{2I} $$ with the total angular momentum $\mathbf{J}$, the moment of inertia $I$. The rotational energies are the corresponding eigenvalues $$ E_\text{rot}/h = BJ(J+1) $$ where $h$ is the Planck constant, $B=\frac{h}{8 \pi^2 I}$ the rotational constant, and $J$ the total angular momentum quantum number. $I$ and $B$ are unique for each molecule and are thus often compared to a molecular fingerprint that can uniquely identify a molecule. The energy levels are shown on the left hand side of the figure below. 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 corresponding 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 below.

The calculated rotational energies (left hand side) and rotational transitions (right hand side) of a linear molecule (OCS at 300 K) up to $J=7$. The energies and frequencies are given in multiples of the rotational constant $B$. In the energy diagram, the arrows indicate transitions and are connected by gray dotted lines with the respective transitions in the spectrum.

When adding centrifugal distortion for a non-rigid molecule, the energy expression expands 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{...} $$ decrease in comparison with the rigid rotor.

The main results are that the rotational spectra of linear molecules have almost equidistant patterns (since $D$ is a few magnitudes smaller than $B$), and that neglecting higher centrifugal distortion constants results in transition frequencies being predicted to slightly higher transition frequencies.

Approximation of the Rotational Constant $B$

Before we examine 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}$. These can be approximated using 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}$).

First, 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 exactly true but is a good first approximation).

Solution

For a two-atomic molecule with atom masses $m_1$ and $m_2$ as well as a distance $r$ between the two atoms, 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 $B$ is given in units of a frequency as $$ B = \frac{h}{8 \pi^2 I} $$ Thus, the radius $r$ for such a two atomic molecule is $$ r = \sqrt{\frac{h}{8 \pi^2 \mu B}} $$

For our case, the masses are $m_\text{C} = 12$, $m_\text{0} = 16$, and $m_\text{S} = 32$. The resulting bond distances for CO and CS are $r_\text{CO} = 1.1308 \cdot 10^{-10} \text{ m}$ and $r_\text{CS} = 1.5375 \cdot 10^{-10} \text{ m}$. We can then estimate the rotational constant of OCS via its moment of inertia: $B=6318.8379 \text{ MHz}$.

The Python implementation of this 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 synthetic spectrum. This has a few educational advantages. Download the spectrum here. It is in *.csv format, with frequencies in MHz in the first column and 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 using either drag-and-drop or the menu (Files > Add Files).

Navigating the Spectrum

First, we want to see an overview of the spectrum. Set the number of plots to one (press Ctrl + N, then enter 1 in the dialog), the center of the plot to 100 GHz (press Ctrl + G, then enter 100000 in the dialog), and the width of the plot to 200 GHz (press Ctrl + W, then enter 2E5 in the dialog). The plot should now look as shown in the figure below.

Overview of the synthetic OCS spectrum in LLWP.

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

Try them out and zoom into some individual peeks.

Loomis-Wood Plots

To see how well our estimated $B$ value compares to the experimental spectrum, we will now switch to a Loomis-Wood representation of the spectrum. Loomis-Wood plots show excerpts of the same width at the predicted positions of the transitions. These excerpts are aligned vertically. If the predictions are perfect, the experimental lines will form a vertical line. Any deviations from the predicted positions can easily be spotted as the experimental lines will form a trend.

Here, we will look at five transitions: $J = 1 \leftarrow 0$ to $J=5 \leftarrow 4$. To do so, follow these steps:

Change the number of plots to 5 (press Ctrl + N, then enter 5 in the dialog)
Change the relative offset of the plots to 0 (press Ctrl + G, then enter 0 in the dialog)
Change the width of each plot to 10 GHz (press 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 (press 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 figure below.

The first Loomis-Wood plot you should see after following this tutorial so far.

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 you entered into the Reference Series window.

The experimental spectrum shifts towards lower frequencies (to the left) as the $J$ values increase (from bottom to top). This clearly indicates that our estimated $B$ value is too high.

We will improve the value of $B$ in the following.

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 width back to 10 GHz (Ctrl + W then enter 10E3 in the dialog). As you can see, the new $B$ value models the series much better and the series almost form a vertical trend. Lets zoom into the Loomis-Wood plot to see that the description of the series is still 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 pronounced, switch to 10 plots and a width of 15 MHz.

The deviations are a result of the centrifugal distortion which we have omitted so-far in our very simple model. As described in the theoretical section, centrifugal distortion parameters are 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 transitions to the correct quantum numbers and then fit the Hamiltonian to these assignments. Based on this Hamiltonian, we would then create new predictions. In the next section, we will learn how to perform these tasks using SPFIT and SPCAT for a simple linear molecule.

Prerequisites

First Steps with SPFIT/SPCAT