Isotopologues

In the previous section, we successfully assigned and modeled the rotational spectrum of the main isotopologue of OCS. Now, we will examine which other isotopologues are present in the spectrum. Isotopologues are very interesting because they are excellent tracers of isotope ratios (e.g., for astronomers) and can be used to determine a molecule's structure. We will start with a short exercise. First, research the natural abundances and masses of the different oxygen, carbon, and sulfur isotopes. Then, calculate the natural abundances of the resulting isotopologues and estimate their rotational constants $B$ in the same way that we did for the main isotopologue previously.

Solution
We copy the values for the $r_{CO}$ and $r_{CS}$ distances from the previous exercise and research the masses and abundances of the most prevalent isotopes. Then, we calculate the abundances of the isotopologues as the product of the abundances of the respective isotopes. For the $B$ rotational constant, we use the same equation as in the previous exercise, but with the masses of the isotopologues. See the following code for a possible Python implementation. 
			
import numpy as np 

r_CO = 1.1308051672090703e-10
r_CS = 1.5375286554975936e-10
m0 = 1.66054e-27
h = 6.626e-34

O_isotopes = [[16, 0.998,   15.994914619257],
				[17, 0.00038, 16.999131755953],
				[18, 0.00205, 17.999159612136]]

C_isotopes = [[12, 0.989, 12],
				[13, 0.0106, 13.003354835336]]

S_isotopes = [[32, 0.948, 31.9720711744],
				[33, 0.0076, 32.9714589099],
				[34, 0.0437, 33.96786701],
				[36, 0.0002, 35.96708070]]


Bs = []
for nO, pO, mO in O_isotopes:
	for nC, pC, mC in C_isotopes:
			for nS, pS, mS in S_isotopes:
				masses = np.array((mO, mC, mS))
				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)
				abundance = pO * pC * pS * 100
				Bs.append(((nO, nC, nS), B/1E6, abundance))


Bs = sorted(Bs, key=lambda x: x[-1], reverse=True)
output = ['| mO | mC | mS | B [MHz] | Abundance|', '|----|----|----|---------|----------|']
for (nO, nC, nS), B, percentage in Bs:
	output.append(f'| {nO} | {nC} | {nS} | {B:.2f} | {percentage:7.4f}% |')

print("\n".join(output))

The results are summarized in the following table:

mO mC mS B [MHz] Abundance
16 12 32 6322.22 93.5697%
16 12 34 6167.36 4.3133%
16 13 32 6301.34 1.0029%
16 12 33 6242.45 0.7501%
18 12 32 5930.97 0.1922%
16 13 34 6144.89 0.0462%
17 12 32 6116.73 0.0356%
16 12 36 6028.68 0.0197%
18 12 34 5780.04 0.0089%
16 13 33 6220.75 0.0080%
18 13 32 5916.16 0.0021%
17 12 34 5963.94 0.0016%
18 12 33 5853.22 0.0015%
17 13 32 6099.16 0.0004%
17 12 33 6038.03 0.0003%
16 13 36 6004.73 0.0002%
18 13 34 5763.87 0.0001%
18 12 36 5644.86 0.0000%
17 13 34 5944.91 0.0000%
18 13 33 5837.72 0.0000%
17 12 36 5827.10 0.0000%
17 13 33 6019.71 0.0000%
18 13 36 5627.42 0.0000%
17 13 36 5806.70 0.0000%

Obviously, we begin our search for isotopologue patterns with the most abundant isotopologue OC$^{34}$S, which is 4%. Its estimated rotational constant is $B=6167.36 \text{ MHz}$ but we have some additional information. From our analysis of the main isotopologue, we know that the estimated rotational constant was slightly too large. Since the rotational constant depends only on the geometry and the masses, our geometry was probably slightly off (which is not very surprising given our estimations). We can calculate a correction factor from the experimental value and the initial estimated value for $B$. This should yield $B_\text{exp} / B_\text{est} \approx 0.962$. Therefore, we should multiply all estimated $B$ values by this correction factor.

Solution
mO mC mS B [MHz] Abundance
16 12 32 6081.98 93.5697%
16 12 34 5933.00 4.3133%
16 13 32 6061.89 1.0029%
16 12 33 6005.24 0.7501%
18 12 32 5705.59 0.1922%
16 13 34 5911.38 0.0462%
17 12 32 5884.29 0.0356%
16 12 36 5799.59 0.0197%
18 12 34 5560.40 0.0089%
16 13 33 5984.37 0.0080%
18 13 32 5691.34 0.0021%
17 12 34 5737.31 0.0016%
18 12 33 5630.80 0.0015%
17 13 32 5867.40 0.0004%
17 12 33 5808.58 0.0003%
16 13 36 5776.55 0.0002%
18 13 34 5544.84 0.0001%
18 12 36 5430.35 0.0000%
17 13 34 5719.00 0.0000%
18 13 33 5615.89 0.0000%
17 12 36 5605.67 0.0000%
17 13 33 5790.96 0.0000%
18 13 36 5413.58 0.0000%
17 13 36 5586.05 0.0000%

Next, we will create a model for the OC$^{34}$S isotopologue. To do so, we can copy the main isotopologue *.par file and only change the $B$ value. It is often better to fix the higher-order constants to the respective main isotopologue values than to set them to zero. You can find the respective *.var file here, in case you want to compare your version. For now, we will use the same *.int file and neglect any small changes in the dipole moment. Create predictions for the transitions by running SPCAT, then load them together with the spectrum into LLWP. You should immediately see a pattern of experimental lines close to the predicted transitions.

The scaled predictions for the OC$^{34}$S isotopologue reproduce the spectrum of OC$^{34}$S already very well. Remaining deviations are easy to follow in the Loomis-Wood plot which makes the assignments efficient and reliable.

As an exercise, create a second set of predictions without the higher-order parameters and compare the two Loomis-Wood plots with each other. You should see that fixing the constants to the main isotopologue value results in more accurate predictions.

Same Loomis-Wood plot as for the previous figure, but without the higher-order constants fixed to the main isotopologue values. The deviations are much more pronounced here (note the different widths of the two plots: 20 MHz for the previous plot and 80 MHz for this one).

Assign-Fit-Predict

Now repeat the same steps as for the main isotopologue:

  1. Assign the correct quantum numbers to the experimental lines in LLWP
  2. Save the assignments to a *.lin file
  3. Create a *.par file from the *.var file (update the uncertainties, and the number of lines)
  4. Run SPFIT
  5. Check the relative uncertainties of the parameters (are all parameters defined, if not, remove them, repeat)
  6. Run SPCAT
  7. Check the residuals

Similarly to the main isotopologue, you can determine sensible uncertainties for your new assignments and adjust the uncertainties in the *.lin file accordingly. Then, search for OC$^{34}$S in the CDMS and add the available literature data to your model. Check if you should add any higher-order parameters. Finally, compare your model to the one found in the CDMS.

No Model - No Problem

After successfully modeling OCS and OC$^{34}$S, we will continue searching for more isotopologues in the spectrum. This time, however, we will forego initial predictions and search for patterns in the Loomis-Wood plots.

Load the OCS spectrum and the *.cat and *.int files of the two isotopologues that have already been modelled. Loading the *.lin files will mark all already assigned lines with a blue star and makes it easy to exclude them from our search. We can then use the predictions for the main isotopologue and search for patterns around it. However, there is one problem. We have loaded two *.cat and *.int files with the same quantum numbers. LLWP uses the first matching prediction for the reference series, which, in this case, are the predictions for OC$^{34}$S. To limit the selection to a single file, go to the Reference Series window and click the All button behind the File: label. In the dialog that opens, choose the file with the main isotopologue . This is very useful when working with many single-state fits, as we are doing here.

Loomis-Wood plot after loading the *.cat and *.lin files of the two isotopologues that have already been analyzed. Red sticks indicate the reference series, the main isotopologue data, while purple sticks indicates any other predicted transitions. Blue stars indicate an assignment.

It is difficult to see any additional isotopologues in the resulting Loomis-Wood plot because each row is scaled to its maximum intensity. Thus, any less intense lines are difficult to spot. We can change this by opening the Scaling Window via the menu (Plot > Scaling Window). In there, we choose the Custom option for the y-axis. New input fields will appear, specifying the minimum and maximum values of the y-axis (excluding a user-defined margin) and the scaling factor between the predicted and experimental intensities. Choose 6.5e2 for the scaling, and enter 0 and 0.001 for the min and max values, respectively. Now, two additional series should be visible in the Loomis-Wood plot.

Two so-far unassigned series become visible in the Loomis-Wood plot when adjusting the y-axis scale.
Awesome Tip
You can align each of the two unassigned series to the center of the plot by using a dynamic offset. Open the offset dialog (Ctrl + G) and enter the following into the Expression field: 
			
-2 * 19.55 * (qnl1 + 1)
LLWP labels the quantum numbers internally as qnu1, qnu2, qnu3, ... for the upper quantum numbers and qnl1, qnl2, qnl3, ... for the lower quantum numbers. Therefore, the expression can be translated to $$ -2 \cdot \Delta B \cdot (J+1) $$ which effectively decreases the $B$ rotational constant by 19.55 MHz.
For the second unassigned series, you can set the $\Delta B$ to about -76.6 MHz. Based on these observations and the previously calculated rotational constants $B$, you can already identify the correct isotopologues for these two series.

Start to assign one of the two series, save the assignments to a new *.lin file, and fit a model to these assignments. Repeat these steps for the second series. Once you have finalized the models for the two series, compare the resulting $B$ values with the estimated $B$ values. Which two isotopologues did you assign?

Solution
You have assigned OC$^{33}$S and O$^{13}$CS. This can be deduced from comparing the predicted $B$ values with the experimental $B$ values or by using physical arguments. The OC$^{33}$S series lies right in the middle of the series for OC$^{32}$S and OC$^{34}$S. For O$^{13}$CS, the change in mass occurs in the carbon, atom which is very close to the center of mass. Therefore, only a small change in the moment of inertia and the rotational constant is expected.

This exercise demonstrates that similar fingerprints (of isotopologues or vibrationally excited states) can be found by searching for similar patterns around already identified patterns in Loomis-Wood plots. Conversely, if you spot one of these patterns in a Loomis-Wood plots, you now know that it probabaly belongs to an isotopologue or a vibrationally excited state.

A good way to practice what you learned in this chapter is to identify isotopologues in the experimental spectrum from the previous chapter.

Linear Molecules in SPFIT/SPCAT
Vibrational States