A Large-scale Approach to Modeling Molecular Biosignatures: The Diatomics

When modeling the rovibrational spectrum, we considered two key parameters: the band origin and the rotational constants. The band origin crucially determines the position of its corresponding rotational lines and thus has a strong influence on the overall shape of the rovibrational band. In our model, we have assumed this aspect to be the most important.

The rotational constants have an effect on the spacing of the transitions on each of the rotational branches. As a secondary effect, anharmonicity causes the rotational constant to depend on the vibrational state. Indeed, the bond lengthening in excited vibrational states allows a Q-branch progression (a set of purely vibrational transitions with no rotational transitions; see Equation (28)), since the rotational constants of both starting and final vibrational states are different. An identical rotational constant for all vibrational level (such as predicted by a simple harmonic oscillator approximation, for example) incorrectly suggests that vibrations and rotations are independent. This leads to the Q--branch transitions bunching up at a single position, rather than displaying a typical progression. Therefore, for an accurate approximation for the entirety of the spectrum, the rotational constant will need to vary with vibrational level (Pawłowski et al. 2002).

The TOSH theory is a modern approach to treating the nuclear wave function that does predict correctly most band origins (Lin et al. 2008). TOSH is effectively a simplification of second-order vibrational perturbation theory (VPT2), and applies limited-order perturbation theory to the nuclear Schrödinger equation. By doing this, it avoids any degeneracy issue that may arise from standard second-order vibrational perturbation theory (Willetts et al. 1990; Vázquez & Stanton 2006), for example.

Fundamentally TOSH works by introducing a shifting parameter (σ) which is used to shift the harmonic basis functions from their equilibrium position. This shift is optimized for the vibrational transition energy expansion, specifically for the fundamental vibrational transition, and as such does an excellent job at correcting the band origins (this is shown by our results in Section 3 and Table 2). An expression for the magnitude of the shift is described in Equations (2)–(10) in Section 2.2. We have discovered, however, that this shift can be applied to the equilibrium distance to recover the anharmonic effects of the first vibrationally excited level (see Section 2.6).

Despite its approximate nature, TOSH is roughly a factor of two cheaper than VPT2, with only some loss in accuracy as a compromise (Lin et al. 2008). This approach also has the potential to be of reasonable accuracy for a large ensemble of molecules (Hanson-Heine 2019) and is straightforward to extend to larger systems (Lin et al. 2008).

We will now briefly discuss the necessary theory required to understand the process to produce the Prometheus spectra.

For a diatomic molecule, the Hamiltonian with up to fourth-order terms (quartic potential) and using mass-weighted displacement coordinate, $Q=\sqrt{\mu }x$, can be expressed as

$$\begin{eqnarray}&&\hat{H}=-\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}}{\partial {Q}^{2}}+\displaystyle \frac{1}{2!}{\eta }_{{ii}}{Q}^{2}+\displaystyle \frac{1}{3!}{\eta }_{{iii}}{Q}^{3}+\displaystyle \frac{1}{4!}{\eta }_{{iiii}}{Q}^{4}.\end{eqnarray} \tag{ 1 }$$

Now imagine a shift, σ, along the coordinate Q. If the center of the wave function is shifted by σ, the shape will remain the same but the anharmonic correction can be incorporated into the wave function.

The shifted wave functions for TOSH are now different to the harmonic wave functions, and are described by

$$\begin{eqnarray}&&{\psi }_{n}={\left(\displaystyle \frac{{\omega }^{1/2}}{{\pi }^{1/2}{2}^{n}n!}\right)}^{1/2}{e}^{-\tfrac{\omega }{2}{(Q-\sigma )}^{2}}{H}_{n}\left[(Q-\sigma ){\omega }^{1/2}\right],\end{eqnarray} \tag{ 2 }$$

where H_n (x) naturally refers to the Hermite polynomial, and ω describes the harmonic frequency (defined using TOSH constants in Equation (13); see below).

The energy of this ground vibrational state is

$$\begin{eqnarray}&&{E}_{0}=\langle {\psi }_{0}| \hat{H}| {\psi }_{0}\rangle ,\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}\begin{array}{rcl}{E}_{0} & = & \displaystyle \frac{1}{4}\omega +\displaystyle \frac{1}{2!}{\eta }_{{ii}}\left[\displaystyle \frac{1}{2\omega }+{\sigma }^{2}\right]+\displaystyle \frac{1}{3!}{\eta }_{{iii}}\left[\displaystyle \frac{3\sigma }{2\omega }+{\sigma }^{3}\right]\\ & & +\ \displaystyle \frac{1}{4!}{\eta }_{{iiii}}\left[\displaystyle \frac{3}{4{\omega }^{2}}+\displaystyle \frac{3{\sigma }^{2}}{\omega }+{\sigma }^{4}\right].\end{array}\end{eqnarray} \tag{ 4 }$$

The equation for the energy of the first excited vibrational level is

$$\begin{eqnarray}&&{E}_{1}=\langle {\psi }_{1}| \hat{H}| {\psi }_{1}\rangle ,\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}\begin{array}{rcl}{E}_{1} & = & \displaystyle \frac{3}{4}\omega +\displaystyle \frac{1}{2!}\left[\displaystyle \frac{3}{2\omega }+{\sigma }^{2}\right]{\eta }_{{ii}}+\displaystyle \frac{1}{3!}{\eta }_{{iii}}\left[\displaystyle \frac{9\sigma }{2\omega }+{\sigma }^{3}\right]\\ & & +\ \displaystyle \frac{1}{4!}{\eta }_{{iiii}}\left[\displaystyle \frac{15}{4{\omega }^{2}}+\displaystyle \frac{9{\sigma }^{2}}{\omega }+{\sigma }^{4}\right].\end{array}\end{eqnarray} \tag{ 6 }$$

Therefore the energy difference between the first vibrational state and the ground state is

$$\begin{eqnarray}&&{\rm{\Delta }}{E}^{\mathrm{TOSH}}=\langle {\psi }_{1}^{* }| \hat{H}| {\psi }_{1}\rangle -\langle {\psi }_{0}^{* }| \hat{H}| {\psi }_{0}\rangle ,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&=\ \omega +\displaystyle \frac{{\eta }_{{iiii}}}{8{\omega }^{2}}+\displaystyle \frac{{\eta }_{{iii}}\sigma }{2\omega }+\displaystyle \frac{{\eta }_{{iiii}}{\sigma }^{2}}{4\omega }.\end{eqnarray} \tag{ 8 }$$

Which can now be compared with the energy from the unshifted wave function obtained through second-order perturbation theory, VPT2:

$$\begin{eqnarray}&&{\rm{\Delta }}{E}^{\mathrm{VPT}2}=\omega +\displaystyle \frac{{\eta }_{{iiii}}}{8{\omega }^{2}}-\displaystyle \frac{5{\eta }_{{iii}}^{2}}{24{\omega }^{4}}-\displaystyle \frac{{\eta }_{{iiii}}^{2}}{32{\omega }^{4}}.\end{eqnarray} \tag{ 9 }$$

Within the TOSH theory derivations, σ is assumed to be small, and therefore the σ² term can be neglected. By comparing the coefficient of η_iii in both the TOSH (8) and VPT2 (9) expressions, a suitable value for the shift can be obtained:

$$\begin{eqnarray}&&\sigma =-\displaystyle \frac{5}{12}\displaystyle \frac{{\eta }_{{iii}}}{{\omega }^{3}}.\end{eqnarray} \tag{ 10 }$$

Due to its derivation, this value of the shift parameter is only optimal for the 0 → 1 transition. It is also worth noting that the original TOSH paper of Lin et al. (2008) does not match the η_iiii terms for the VPT2 energy expression, but remedied this by stating this term is often neglected, in any case. Many intermediate processes for the derivations above have been assumed and therefore omitted. For a full derivation, please see the Appendix, Section A.1.

One lesser-considered aspect of having a displaced vibrational wave function (mimicking what happens for the exact wave function) is the usage the shift parameter in order to include anharmonic effects in the computation of rotational constants. We explore this possibility in Section 2.6.

As shown in Equations (8) and (10), the TOSH approach (and our Prometheus implementation; see Section 2.8) only requires a selection of anharmonic constants for the chosen diatomic. These constants can be obtained either through finite-difference methods or by computing the derivatives of the potential energy curve (PEC) directly, if a functional form is available. In this study, we do not focus on producing PECs, but we use existing ones to derive the required anharmonic constants using a quartic fitting procedure described below.

Prometheus uses the numerical points that make up the PEC to perform a quartic fit (see Equation (11)) and extract the anharmonic constants. We use this methodology for two reasons. Historically quartic force fields have been used as a generic polynomial form to model PECs. In particular, it is used as a tool for analyzing and producing rovibrational spectra for molecules of interest to astrophysical observation (Fortenberry & Lee 2019). Second, the TOSH theory requires only the second-, third- and fourth-order derivatives, without any need for higher-order derivatives.

The quality of the PEC is paramount and will affect the accuracy of the TOSH constants derived from our quartic fit. The code performs a local fit, centered around the equilibrium bond length, rather than a global fit, which would include the potential wall. The local fit is determined via an "inclusion range" determined by starting with the equilibrium bond length and including points at differing percentages of the equilibrium bond length. The ranges tested within the code span from ±20% to ±7% of the equilibrium bond length, in intervals of 1%.

The fit accuracy has been tested for each range using a simple chi-square test. The range with the lowest chi-square value, or the first range with a sufficiently low chi-square value (arbitrarily defined as 1 × 10⁻⁷) was selected. This was to ensure the range providing the best fit (while still considering the point density of the remaining curve) to the reference data is selected and used to ultimately produce a rovibrational spectrum. This method has been chosen to ensure a consistent selection of inclusion ranges for the PEC of each molecule. It also removes any bias we could have introduced into the analysis by selecting the ranges ourselves.

The quartic fit applied to the PEC is described by the equation

$$\begin{eqnarray}&&V={E}_{0}+\displaystyle \frac{1}{2!}{\zeta }_{{ii}}{X}^{2}+\displaystyle \frac{1}{3!}{\zeta }_{{iii}}{X}^{3}+\displaystyle \frac{1}{4!}{\zeta }_{{iiii}}{X}^{4},\end{eqnarray} \tag{ 11 }$$

where X = r − r_e is a displacement coordinate, with r_e describing the equilibrium bond length. Note that E₀ has been included to help Prometheus with the fitting by centering around 0, and has little importance beyond this.

The fit only returns ζ constants, ζ_ii , ζ_iii , and ζ_iiii, which need to be massed weighted for them to become the TOSH constants, as described in the Hamiltonian in Equation (1). This was done in the following manner:

$$\begin{eqnarray}&&{\eta }_{{ii}}=\displaystyle \frac{{\zeta }_{{ii}}}{\mu },\quad {\eta }_{{iii}}=\displaystyle \frac{{\zeta }_{{iii}}}{{\mu }^{3/2}},\quad {\eta }_{{iiii}}=\displaystyle \frac{{\zeta }_{{iiii}}}{{\mu }^{2}}.\end{eqnarray} \tag{ 12 }$$

Note that fits of order higher than quartic lead to sizable errors and inaccurate results. Additionally, we used spline interpolation to ensure the data spanned the entire range specified.

In this study, we obtain the anharmonicity constants by computing the derivatives of a quartic fit to a given PEC. As previously alluded to, these constants are the only data required for Prometheus to operate. The technique used to derive these constants (finite difference, derivatives of PEC, etc.) are inconsequential, but require a suitable description of the interatomic potential around the equilibrium bond length in order to produce reliable spectra. In order to avoid issues originating from a PEC of subspectroscopic quality, we choose to use only data that have been validated through spectroscopic comparisons.

Indeed, while a number of PECs for diatomics are available in the literature, some of the molecules in our selected set have a particularly challenging electronic structure (CO, N₂ and O₂). This in turn affects the quality of the anharmonicity constants, which then impacts the accuracy of the rovibrational lines. As our study focuses on the quality of the rovibrational approach, rather than the quality of the electronic structure approach used for the potential, we chose well-reported modern PECs, obtained through either Rydberg–Klein–Rees (RKR) inversion (see Rees 1947 for further details) or ab initio calculations. All the potentials used within this study are for the ground state of their respective molecule.

The PECs used are those from Schwartz & Le Roy (1987) (the mass-independent clamped nuclei, V_BO, potential with adiabatic correction, ΔV_ad) for the ${}^{1}{{\rm{\Sigma }}}_{g}^{+}$ state of H₂. The PEC for ¹Σ⁺ for CO is from Meshkov et al. (2018) (specifically, the CO_X_Func_edited Born–Oppenheimer (UBO), potential located within the supplementary material). The PEC for the ${}^{3}{{\rm{\Sigma }}}_{g}^{-}$ state of O₂ is from Yu et al. (2014) and the potential for the ${}^{1}{{\rm{\Sigma }}}_{g}^{+}$ state of N₂ is from Edwards et al. (1993).

H₂ is an ab initio potential, CO a semiempirical potential, while N₂ and O₂ are spectroscopic RKR curves. This has been done in part to highlight the flexibility of our model, which can accept any type of PEC, rather than being bound to a single format. Naturally, if the quality of the curve is low, the rovibrational results will reflect this fact, but will still remain a better approximation than a pure harmonic approach.

The necessary anharmonicity or TOSH (quartic) constants, obtained using optimal quartic fitting ranges, for each molecule in this study are given in Table 1.

Table 1. Anharmonic Constants Computed by Prometheus (TOSH) for Each Molecule

Constant	H2	O2	N2	CO
ηii (${E}_{h}{a}_{0}^{-2}{m}_{e}^{-1}$)	4.040 × 10−4	5.183 × 10−5	1.156 × 10−4	9.803 × 10−5
ηiii (${E}_{h}{a}_{0}^{-3}{m}_{e}^{-3/2}$)	− 4.971 × 10−5	− 1.718 × 10−6	− 4.048 × 10−6	− 3.430 × 10−6
ηiiii (${E}_{h}{a}_{0}^{-4}{m}_{e}^{-2}$)	5.317 × 10−6	4.789 × 10−8	1.069 × 10−7	8.928 × 10−8

Note. Atomic units are used throughout, where E_h , a₀, and m_e represent Hartrees, Bohr radius, and electron rest mass, respectively. The constants are derived from a quartic fit to the numerical PECs (see Section 2.3 for details).

Download table as:  ASCIITypeset image

The TOSH framework allows us to determine key constants, such as the harmonic frequency, ω, the anharmonic fundamental transition, ΔE^TOSH or ν_0→1, and the coordinate shift, σ.

The value of ω is calculated using

$$\begin{eqnarray}&&\omega ={\left({\eta }_{{ii}}\right)}^{1/2}.\end{eqnarray} \tag{ 13 }$$

The value of ΔE^TOSH is obtained from Equation (8), and the shift, σ, is given by Equation (10). The associated errors for each spectroscopic constant is also calculated to allow an assessment of the error on the positions of transitions in the spectra produced by Prometheus (see Section 2.8).

The state-specific rotational constant, B_v , for a diatomic molecule is defined as

$$\begin{eqnarray}&&{B}_{v}=\displaystyle \frac{{h}^{2}}{8{\pi }^{2}c\mu {r}_{v}^{2}},\end{eqnarray} \tag{ 14 }$$

where v refers to the vibrational level, μ is the reduced mass, and r_v is a bond length that depends on the vibrational energy level considered. We can see that the bond constant has a reciprocal squared relationship to the bond length, meaning as length increases the constant will decrease. This is one of the fundamental equations used throughout Prometheus to ultimately create the rovibrational spectra.

In the harmonic approximation, r_v is independent of the vibrational energy level and thus r_v = r_e , the equilibrium bond length. Consequently, in that approximation B_v = B_e , a fixed rotational constant obtained from r_e . More generally, the value of the bond length for the rotational constant for each vibrational state can be obtained from the expectation value of the position, r_v = 〈ψ_n ∣r∣ψ_n 〉. Using this expression for the harmonic model leads to the same conclusions as earlier: r_v = r_e , since the harmonic wave function is symmetric around r_e .

If we now consider the TOSH model, the shifted TOSH position, x, can be described using the harmonic position, u, and the shift constant, σ:

$$\begin{eqnarray}&&x=u+\sigma .\end{eqnarray} \tag{ 15 }$$

As shown earlier, we can take the expectation value of the position, but this time implementing TOSH:

$$\begin{eqnarray}&&\langle {\psi }_{n}^{T}| x| {\psi }_{n}^{T}\rangle =\langle {\psi }_{n}^{H}(u)| u+\sigma | {\psi }_{n}^{H}(u)\rangle ,\end{eqnarray} \tag{ 16 }$$

where ${\psi }_{n}^{T}$ and ${\psi }_{n}^{H}$ represent the TOSH (Equation (2)) and harmonic wave functions, respectively:

$$\begin{eqnarray}&&\langle {\psi }_{n}^{T}| x| {\psi }_{n}^{T}\rangle =\langle {\psi }_{n}^{H}(u)| u| {\psi }_{n}^{H}(u)\rangle +\langle {\psi }_{n}^{H}(u)| \sigma | {\psi }_{n}^{H}(u)\rangle .\end{eqnarray} \tag{ 17 }$$

The first expectation value is the position expectation position for a harmonic wave function (i.e., zero) and thus we can see that the expectation of the position for TOSH is simply

$$\begin{eqnarray}&&\langle {\psi }_{n}^{T}| x| {\psi }_{n}^{T}\rangle =\sigma .\end{eqnarray} \tag{ 18 }$$

The full derivation, for both required levels, can be found in the Appendix, Section A.2.

Thus the TOSH model leads to r_v = r_e + σ. As is the case for the harmonic model, the TOSH approach only produces a single static value for all rotational constants (see Section A.2).

We cannot accurately reproduce rotational constants that exhibit rotation–vibration interaction, using solely either a TOSH or a harmonic model: a hybrid theory is required. The approach we use in Prometheus approximates the ground-state rotational constant using the harmonic value. The first excited vibrational level uses the modified TOSH bond length derived above, where sigma is combined with the equilibrium bond length. The vibrationally averaged bond distances are thus expressed as follows:

$$\begin{eqnarray}&&{r}_{0}={r}_{e},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{r}_{1}={r}_{e}+\sigma ,\end{eqnarray} \tag{ 20 }$$

where σ is not mass weighted to harmonise units.

The harmonic approach effectively costs nothing (using r₀ = r_e ), so can be exploited to provide the ground-state rotational constant with little effort. Our approach does not necessarily provide a strong quantitative agreement but is qualitatively correct: bond length gets larger with increased vibration, hence the rotational constant changes (gets smaller).

In order to generate a rovibrational line list and corresponding spectrum, we compute the maximum allowed quantum rotational number, J, for a given input temperature. In the present study, the temperature was set to 300 K to match the other data sets. The spectra from each of the experimental databases are all set at 300 K, also.

Rotation symmetry also influences the spectrum the molecule will produce, and therefore conditional arguments within the intensity calculations have been created to account for this. If a diatomic molecule has a center of symmetry, hence meaning it is homonuclear, it is Raman active. This occurs instead of IR activity due to having no dipole moment; instead stretching and contraction of the bond leads to changes in the polarizability. This type of molecular response implies different transition selection rules leading to the O, Q, and S branches, rather than the typical P and R branches in infrared.

Nuclear spin also has an effect on the spectra produced. This will be discussed further in Section 3. The code is automated to account for additional effects by using the initial inputs of the mass of a molecule.

We use a Boltzmann distribution to obtain the relative intensities for each transition. This statistical law states that, for a system of N total molecules, only a fraction, $\tfrac{{N}_{J}}{N}$, will occupy particular energy level, E_J , a fraction. This can be written as

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{J}}{N}=\displaystyle \frac{{g}_{J}{e}^{\left(\tfrac{-{E}_{J}}{{kT}}\right)}}{f},\end{eqnarray} \tag{ 21 }$$

where the degeneracy, g_J , the energy of a given level, E_J , and the partition function, f, are described as

$$\begin{eqnarray}&&{g}_{J}=(2J+1),\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{E}_{J}={B}_{0}J(J+1),\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&f=\sum {g}_{J}{e}^{\left(\tfrac{-{E}_{J}}{{kT}}\right)}.\end{eqnarray} \tag{ 24 }$$

Hence, the full equation now becomes

$$\begin{eqnarray}&&\displaystyle \frac{{N}_{J}}{N}=\displaystyle \frac{(2J+1){e}^{\left(\tfrac{-{BJ}(J+1)}{{kT}}\right)}}{f}.\end{eqnarray} \tag{ 25 }$$

From here the code calculates the transition positions using the appropriate equations depending on the rotational branches considered. The full derivations for these equations have been included in the Appendix, Section A.3. They are summarized as follows:

$$\begin{eqnarray}&&{\nu }_{O}={\nu }_{0\to 1}+{B}_{1}J^{\prime} (J^{\prime} +1)-{B}_{0}(J^{\prime} +2)(J^{\prime} +3),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\nu }_{P}={\nu }_{0\to 1}-{({B}_{1}+{B}_{0})(J^{\prime} +1)+({B}_{1}-{B}_{0})(J^{\prime} +1)}^{2},\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{\nu }_{Q}={\nu }_{0\to 1}-J^{\prime} (J^{\prime} +1)({B}_{1}-{B}_{0}),\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\nu }_{R}={\nu }_{0\to 1}+({B}_{1}+{B}_{0})(J^{\prime\prime} +1)+({B}_{1}-{B}_{0}){\left(J^{\prime\prime} +1\right)}^{2},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\nu }_{S}={\nu }_{0\to 1}+{B}_{1}(J^{\prime\prime} +2)(J^{\prime\prime} +3)-{B}_{0}J^{\prime\prime} (J^{\prime\prime} +2),\end{eqnarray} \tag{ 30 }$$

where ν_0→1 represents the band origin, and B₀ and B₁ represent the rotational constants for their appropriate vibrational levels. As previously mentioned, for our code we use the solution calculated by TOSH (Equation (8)) as the band origin.

Here is a reminder of the labeling for the line series: ΔJ = − 2 (O branch), ΔJ = − 1 (P branch), ΔJ=0 (Q branch), ΔJ = + 1 (R branch), and ΔJ = + 2 (S branch).

Finally, the code combines the intensity calculations with the transition positions to produce the synthetic spectra.

Prometheus is the name given to our code, which is written in Python. It is maintained by the authors and is available on the Milne Centre Github (https://github.com/Milne-Centre/Prometheus). The version of Prometheus used to obtain the results of this paper is archived at 10.5281/zenodo.5494420.

Our code is designed to be an approximate method of simulating spectra, potentially trading off accuracy for simplicity and speed. The main idea behind this code is to offer a new complementary approach to modeling the vast amount of molecules which have astrophysical and astrobiological importance.

Our analysis may have applications beyond the SBP catalog, as we have done ourselves within this work by modeling CO.

Prometheus primarily implements the TOSH theory detailed in Section 2.2 and makes use of the shift parameter, σ, to calculate band origins and provide anharmonic corrections to spectroscopic constants. To our knowledge this is the first time it has been used to provide anharmonic corrections to rotational constants, rather than purely being used for anharmonically corrected vibrational frequency calculations.

Prometheus then produces a stick spectrum (without any line-broadening effects) using the constants it has determined.

The first biosignature molecule is molecular hydrogen. H₂ is the most abundant molecule in the universe by orders of magnitude (Wakelam et al. 2017). It is difficult to observe directly in the interstellar medium due to a lack of permanent electric dipole moment and most transitions are of quadrupolar nature (and extremely weak) (Herzberg 1949; Roueff et al. 2019), or in emission for selected environments (see Habart et al. 2005 for a review). Indirect techniques, such as the amount of dust present using the gas-to-dust ratio (Jo et al. 2017), also allow some degree of detection. See Seager et al. (2020) for further discussion of possible H₂ presence and detection on exoplanets. From a biology standpoint, many microorganisms produce hydrogen as primary product through metabolism under anaerobic conditions (Schlegel 1974).

For H₂, data from HITRAN2016 (Gordon et al. 2017) were available (specifically, the line list comes from Komasa et al. 2011 and Wolniewicz et al. 1998). The spectroscopic constants were obtained from Campargue et al. (2012). Additionally, spectra were available from ExoMol (Roueff et al. 2019).

In Table 2, Campargue et al. (2012) determined the band origin to be 4161.17 cm⁻¹. The harmonic approach gives 4411.45 cm⁻¹, which deviates from the literature by 250.28 cm⁻¹. TOSH, with a result of 4174.71 cm⁻¹, manages to produce a band origin with a difference of only 13.54 cm⁻¹, an order of magnitude less than the harmonic approach.

The experimental ground-state rotational constant, B₀, was determined by Campargue et al. (2012) to be 59.3329 cm⁻¹. The harmonic method approximates B₀ better, with a value of 60.9513 cm⁻¹, than TOSH, which appears to drastically overcorrect, leading to a value of 54.2350 cm⁻¹.

However, since the predicted rotational constants has no vibrational dependence for TOSH or the harmonic model, the reverse is true for B₁. Campargue et al. (2012) determined 56.3732 cm⁻¹, with TOSH and harmonic now differing by 2.1382 cm⁻¹ and 4.5781 cm⁻¹, respectively.

Figures 2 and 3 show the resulting simulated spectra from Prometheus compared to ExoMol, data by Campargue et al. (2012), and the harmonic method.

Zoom In Zoom Out Reset image size

Two major effects need to be considered when the rovibrational spectrum of H₂ is analyzed. First, this homonuclear molecule exhibits a center of symmetry. As a consequence of this, it is Raman active and IR inactive with O, Q and S branches present in its spectra (Campargue et al. 2012). Prometheus will model a molecule in the correct region by considering symmetry of the diatomic using the initial mass inputs.

As also mentioned by Campargue et al. (2012), a molecule composed of nuclei with nonzero nuclear spin (such as H₂) will exhibit spectral lines that show an alternation in intensity. Unfortunately Prometheus does not currently model this. Instead the intensities have the typical Boltzmann distribution regardless of the nuclear spin of the molecule's nuclei components.

A key aspect of Figure 2, is how well all three band origins align. Using spectroscopic constants of Campargue et al. (2012), Prometheus can quite accurately reproduce the ExoMol positions for the transitions. Prometheus TOSH approach however appears to slightly overestimate the band origin by about 20 cm⁻¹. This is likely due to the strong quantum nature of the hydrogen nuclei, and thus this molecule displays strong anharmonic effects that are less well modeled by the TOSH approach. Regardless, Figure 2 shows all three approaches are in relatively good agreement with one another.

Regarding the Q branch, Prometheus displays a larger spread of transitions than ExoMol and Campargue et al. (2012). The Prometheus central transitions potentially would roughly line up provided the band origin was corrected. Once again Campargue et al. (2012) adequately models the ExoMol data. It does falter slightly with the distribution of the branch, as the Campargue et al. (2012) Q branch appears to model a greater spread in the transitions than ExoMol.

The O branch, indicated by the labels on Figures 2 and 3, shows an interesting result. Campargue et al. (2012) appears to decrease slightly in accuracy the higher the transition, whereas Prometheus does the same but the decrease in accuracy is larger. This deterioration is to such a point it becomes difficult to match the ExoMol's transitions to Prometheus.

Zoom In Zoom Out Reset image size

For the S branch (also indicated by the label on Figures 2 and 3), Prometheus seems to initially fare quite well. The Prometheus transition positions appear to slightly increase in difference to the ExoMol data with each higher transition. Possibly if the origins had aligned better Prometheus would be able to continue to quite effectively model this branch at higher rotational transitions. As before, Campargue et al. (2012) gives a satisfactory level of comparison to ExoMol, until the higher transitions are reached at which point it begins to falter. Some difficulty occurs with the comparisons at higher transitions as the relative intensities for the ExoMol data are very low. The ExoMol S branch seems to have a far greater weighting than the O branch, something which is not modeled in any other spectra.

The difficulties in describing the higher rotational transitions can be partly explained by the fact that hydrogen is a light molecule and thus distortion constants are required to accurately model anything beyond J = 0, 1. For example, it was found by Jennings & Brault (1983) that to accurately fit up to J = 5, four rotational constants (B, D, H, and L) were required. This presumably explains Prometheus' inaccuracies for the hydrogen spectra, which are not as pronounced in the other heavier molecules considered within these simulations.

Finally, we compare the Prometheus spectrum with the harmonic model in Figure 3. As discussed earlier, it can be seen that Prometheus slightly overestimates the origin, whereas the harmonic method substantially overestimates it. Due to this, it is difficult to comment on much of the harmonic spectrum, as it barely matches ExoMol.

We can see that despite Prometheus also having some difficulties with modeling H₂, particularly at higher transitions and in the O branch, it still performs better than the harmonic approach. Certainly, this shows anharmonic corrections are required to produce a qualitatively correct result.

The second biosignature is molecular oxygen. Oxygen is the third most abundant element in the universe, yet molecular oxygen is one of the most elusive molecules (Wang et al. 2020). This is often cited to be due to O₂'s high chemical reactivity and lack of electric dipole moment (Luspay-Kuti et al. 2018). Even now, a comprehensive picture of oxygen chemistry in interstellar environments is still missing (Wang et al. 2020; Wakelam et al. 2010; Agundez & Cernicharo 2006).

Oxygen is crucial to our understanding of most life on Earth, with oxygenic photosynthesis being the dominant producing metabolism on our planet (Domagal-Goldman et al. 2014). Concerning exoplanet observations, oxygen is a potential biosignature since it can be produced through photosynthesis processes. However, it is also a possible false positive since it can be formed in larger quantities through abiotic processes such as runaway greenhouse effects (Meadows 2017).

Oxygen's potential to act as a false-positive biosignature is rooted in the diversity of exoplanets. For example, on Earth there are no abiotic processes that would produce it in large abundance (Meadows et al. 2018). It has been shown for a very different star and planetary system that O₂ could be generated in a planetary environment without life (Wordsworth & Pierrehumbert 2014).

The literature spectroscopic constants for O₂ were obtained from Yu et al. (2014), and our values are compared to data from the HITRAN 2016 release (Gamache et al. 1998; Yu et al. 2014; Gordon et al. 2017) only, as no ExoMol data were available. It is worth noting that the HITRAN data contain not Raman but magnetic dipole and electric quadrupole IR transitions, which may lead to some discrepancies between the data sets.

Like H₂, O₂ has a center of symmetry, which causes it to be Raman active but IR inactive. This once again means that the O, Q, and S branches are present in the spectra. In addition, the symmetry also means the effects of nuclear spin will be observed, but these effects differ to that of the other homonuclear diatomic molecules.

Unlike the other diatomic molecules in this study, O₂ has a zero nuclear spin, which has a different effect on the spectra than the alternating intensities effect. Instead, every rotational level with an even value is absent from the spectra. This presents itself as a spectra with seemingly large gaps between the rotational transitions (Yu et al. 2014; see Figure 4).

Zoom In Zoom Out Reset image size

O₂ has a triplet sigma state and exhibits a triplet structure for each of the remaining transition lines, due to the splitting of the rotational levels. This is another effect that Prometheus does not currently model and needs to be taken into consideration when evaluating spectra. However, the weighting of the relative intensities due to triplet splitting have been included into Prometheus, to allow for ease of comparison.

When considering the band origins (see Table 2), TOSH produces a better approximation to the literature than the harmonic method, despite the points mentioned previously. Yu et al. (2014) reports the band origin to be at 1556.39 cm⁻¹, and this is extremely well estimated by TOSH, with a result of 1556.43 cm⁻¹. The harmonic approximation, on the other hand, displays a significant displacement, roughly 24 cm⁻¹. We also note again that the harmonic approximation provides a better approximation of B₀, with TOSH doing the same for B₁.

In Figure 5, we can quite clearly see the shift in band origin has a significant effect on the harmonic approximation's ability to mimic the literature.

Zoom In Zoom Out Reset image size

Prometheus, by using our hybrid methodology, models the Q branch well (in good agreement with Yu et al. 2014). This is something the purely harmonic method (or a purely TOSH approach) repeatedly fails at due to the fixed rotational constants between levels.

The third diatomic biosignature is molecular nitrogen. This is a key species in cosmology and has been observed in various galactic environments (Vangioni et al. 2018), although direct observation of N₂ is often difficult as it lacks strong pure rotational or vibrational lines (Li et al. 2013).

Within our own solar system we find N₂ present in various planetary and satellite atmospheres; it accounts for 3.5% (Oyama et al. 1980) of the Venusian atmosphere, 78.1% (Cox 2002) of the terrestrial atmosphere, 2.8% (Franz et al. 2017) of the Martian atmosphere, and, perhaps most significantly, 98.4% (Strobel & Shemansky 1982) of Titan's atmosphere.

From a biogenic perspective, N₂ is an essential ingredient for the building blocks for life as we know it (Sproß et al. 2018) since it is required, along with carbon and phosphorus, for the formation of nucleic acids and proteins (Lammer et al. 2019).

The spectroscopic constants were obtained from Bendtsen & Rasmussen (2000). Our N₂ results are compared to the harmonic approximation, literature optimal results and the data from the HITRAN 2016 release (Roy et al. 2006; Li & Roy 2007; Gordon et al. 2017). The available ExoMol data was not suitable for our study as it provides the Western et al. (2018) data, which focus on a different region of the spectrum (4500–11,000 cm⁻¹), rather than the ground state.

Finally, N₂, like O₂ and H₂, has a center of symmetry and thus the spectrum is Raman active but IR inactive. As N₂ is composed of nuclei of nonzero nuclear spin, the HITRAN spectrum shows the characteristic alternation in intensities (again, not modeled here by our Prometheus approach).

By looking at the results of Table 2, we see that N₂ rotational constants follow the typical pattern of the previous molecules within this study. Bendtsen & Rasmussen (2000) determined the first excited level rotational constant to be 1.9722 cm⁻¹ and the ground rotational constants to be 1.9896 cm⁻¹. The TOSH approach, with a fixed value of 1.9753 cm⁻¹, overestimates the ground state but provides a good approximation for the first excited level. Whereas the harmonic approximation, with a constant value of 1.9982 cm⁻¹, achieves the opposite.

Zoom In Zoom Out Reset image size

We can see in Table 2 that for N₂, TOSH (2329.85 cm⁻¹) once again provides a superior approximation of the band origin than the harmonic approximation (2359.56 cm⁻¹), when comparing to the literature value (2329.91 cm⁻¹). In Figure 7, we can draw comparisons between HITRAN, Prometheus, and the Bendtsen & Rasmussen (2000) constants used in a Prometheus spectrum. As expected from using the TOSH band origin, Prometheus provides an excellent approximation of both Bendtsen & Rasmussen (2000) and HITRAN.

Zoom In Zoom Out Reset image size

For the lower rotational transitions of both the O and S branches, Prometheus does a satisfactory job of modeling the positions. In most cases it lines up incredibly well with not only Bendtsen & Rasmussen (2000) but also HITRAN data. As was the case for the other molecules, this accuracy is lower for higher J transitions.

Prometheus does an excellent job at approximating the Q branch for N₂ and the relative spread of the transitions of this branch is comparable to the HITRAN data.

For each increasing J transition the position is incrementally overestimated, but not to the point where the Prometheus spectrum is incorrectly matching lower number transitions to higher HITRAN ones. With this is mind, Prometheus is a very good approximation of both Bendtsen & Rasmussen (2000) and HITRAN.

A comparison with the harmonic spectrum is shown in Figure 7. The harmonic band origin is displaced from the literature by roughly 30 wavenumbers (see Table 2). This effectively renders any comparisons of the transition positions to the other spectra pointless.

As was the case for H₂ the relative intensities of N₂ alternate in weighting, an effect that Prometheus does not currently model, once again causing the difficulties in drawing comparisons. The general intensities of the harmonic and the Bendtsen & Rasmussen (2000) spectra have the same problem, which is expected, as it is produced via our Prometheus methodology.

Comparing the relative intensities of Prometheus and harmonic to HITRAN does raise interesting points. The "peak" for the O branch for HITRAN appears to occur 10 to 20 wavenumbers away from Prometheus. Potentially the alternating intensities are artificially causing the "peak" to occur at higher wavenumbers. The S branch does not appear to be as exaggerated. As discussed earlier, this discrepancy is not a major concern as intensities calculations are the not the focus of our present study.

We may conclude that, for N₂, Prometheus is producing a better approximate spectrum than the harmonic method. Like the other molecules seen in previous sections, Prometheus' predictive power is reduced at higher rotational transitions. Presumably the uncorrected harmonic value for ground rotational constant and the TOSH approximate values for the rotational constant of the first excited level is having an effect on the later transitions' positions.

The last molecule considered in this work is CO. Unlike the previous molecules, CO has a permanent dipole moment and therefore readily absorbs, even at the low temperatures found in molecular clouds (Whitworth & Jaffa 2018). More relevantly, CO has already been detected in an exoplanet's atmosphere (HD 189733; de Kok et al. 2013).

Within a biogenic context, CO is often classified as an antibiosignature. For example, for an inhabited planet, it could be difficult for CO to accumulate in the atmosphere as it acts as an energy source for some microbes on Earth (Wang et al. 2016). Anthropologically, it can be formed instead when any organic substance is combusted incompletely and therefore can be an indication of not just life but intelligent life (Horner 2000).

The literature spectroscopic constants were obtained from Coxon & Hajigeorgiou (2004). Our results for CO are compared to the harmonic approximation, optimal results, and data from both the HITRAN 2016 release (Coxon & Hajigeorgiou 2004; Li et al. 2015; Gordon et al. 2017; Devi et al. 2018) and ExoMol's Li2015 line list (Li et al. 2015).

CO does not have a center symmetry and therefore it is IR active and Raman inactive. It is easily modeled by Prometheus to an accuracy akin to the HITRAN data at the lower transitions. Unlike the other biogenic diatomic molecules in this analysis, CO solely demonstrates P and R branches, with no Q-branch transitions.

Coxon & Hajigeorgiou (2004), as shown in Table 2, provided 1.9225 cm⁻¹ for B₀ and 1.9050 cm⁻¹ for B₁. As before, the harmonic approximation for B₀ is better (1.9316 cm⁻¹) than TOSH's 1.9080 cm⁻¹. For the upper rotational constant, B₁, TOSH provides a rotational constant closer to literature.

The band origin of the Prometheus spectrum (using TOSH), is 2143.15 cm⁻¹, and easily lies within a single wavenumber of Coxon & Hajigeorgiou (2004) (at 2143.27 cm⁻¹). Figure 8 shows that, even when compared to the HITRAN/ExoMol spectra, Prometheus' band origin still lies within a few wavenumbers.

Zoom In Zoom Out Reset image size

The spectrum produced via the harmonic method (Figure 9) gives a band origin of 2173.07 cm⁻¹, which has a difference of 29.80 cm⁻¹ to Coxon & Hajigeorgiou (2004). It is also displaced from ExoMol/HITRAN by approximately 30 to 40 wavenumbers. Such a displacement naturally renders the harmonic CO spectrum considerably inaccurate when considering the positions of individual J transitions.

Zoom In Zoom Out Reset image size

In general the distributions of the intensities between HITRAN and Prometheus are similar. On closer inspection, only the fundamental transitions at lower rotational energy levels, J, are modeled with high accuracy.

Indeed, most of the lower-target transitions lie within a couple of wavenumbers of Prometheus transitions, but this starts to falter for the higher values of J. The error for the higher transitions is then exaggerated due to the slight position differences of the band origin.

The relative intensities appear to be in a good agreement, although it is worth mentioning that ExoMol is exhibiting some increased intensities for certain higher-level transitions. This is not shown by either HITRAN2016 or Prometheus data.

As shown in Table 2, the harmonic approximation consistently describes the ground level rotational constant (B₀) well but poorly estimates the upper rotational constant (B₁). This makes realistically modeling any Q branch impossible since it relies on a difference between the ground and excited state rotational constant (see Equation (28) and Section 2.7).

The TOSH-corrected rotational constants, rather than offering the predicted "vibrationally averaged" value between the first and ground states, actually model the first excited state well (see Table 2). Unfortunately, the TOSH rotational constant overestimates the ground-state value, with the harmonic approach typically offering a better approximation. Other studies have also shown (e.g., Table 1 in Endres 1967) that the uncorrected harmonic approximation is a good approximation for the ground state, whereas for the higher levels this is not the case.

Both the ground and first-excited-level rotational results have been included in Table 2 and discussed within the previous results sections. This has been done to highlight that using the current modified theory, where we amalgamate the TOSH and harmonic results, is the best option to match reference spectra. The results used for Prometheus are highlighted in the table for each case.

Typically Prometheus is able to replicate the band origin within a single wavenumber using TOSH, with the exception of H₂. Prometheus has also shown it is capable of estimating the band origins to a greater accuracy than that of the harmonic method, often with order of magnitude improvements.

While TOSH can approximate the upper rotational constant better than the harmonic model approximates the lower one, both still only approximate the literature values.

This has the further effect of causing incorrect spacing of transitions in the branches, which has been shown for the high-J rotational transitions in all spectra produced. We can verify this conjecture by reviewing Equations (26)–(30). Each of these position equations have a quadratic dependence on the quantum number, J, therefore any variation between the rotational constants (such as the differences between the literature and Prometheus) will become more apparent at the higher transitions.

It is also worth noting that we have not used any distortion constants or additional rotation constants in our calculations. We did this to keep the method as simple as possible. The effects of this exclusion are most notable for H₂.

We will go through how to derive the corresponding energy for each vibrational level as referenced by Equations (3) and (5). We will use the following notation from Section 2.6 within the following Appendix sections. As a reminder, the relation between the anharmonic position, x, and the harmonic position, u, can be described as

$$\begin{eqnarray}&&x-\sigma =u.\end{eqnarray} \tag{ A1 }$$

We will now demonstrate deriving the ground-state energy. To begin we must use perturbation theory. As the anharmonic Hamiltonian (Equation (1)) and wave functions are known (Equation (2)) we can express it as

$$\begin{eqnarray}&&{E}_{0}=\langle {\psi }_{0}^{* }| \hat{H}| {\psi }_{0}\rangle ,\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&={\int }_{-\infty }^{\infty }{\psi }_{0}^{* }\left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}+\displaystyle \frac{1}{2!}{\eta }_{{ii}}{x}^{2}+\displaystyle \frac{1}{3!}{\eta }_{{iii}}{x}^{3}+\displaystyle \frac{1}{4!}{\eta }_{{iiii}}{x}^{4}\right){\psi }_{0}{dx}.\end{eqnarray} \tag{ A3 }$$

For ease of comprehension, we will break the following derivations down into four components and solve individually, then sum back together. We will now describe the separated parts and label for ease:

$$\begin{eqnarray}&&A={\int }_{-\infty }^{\infty }{\psi }_{0}^{* }\left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\delta }^{2}}{\delta {x}^{2}}\right){\psi }_{0}\,{dx},\end{eqnarray} \tag{ A4 }$$

$$\begin{eqnarray}&&B={\int }_{-\infty }^{\infty }{\psi }_{0}\left(\displaystyle \frac{1}{2!}{\eta }_{{ii}}{x}^{2}\right){\psi }_{0}\,{dx},\end{eqnarray} \tag{ A5 }$$

$$\begin{eqnarray}&&C={\int }_{-\infty }^{\infty }{\psi }_{0}\left(\displaystyle \frac{1}{3!}{\eta }_{{iii}}{x}^{3}\right){\psi }_{0}\,{dx},\end{eqnarray} \tag{ A6 }$$

$$\begin{eqnarray}&&D={\int }_{-\infty }^{\infty }{\psi }_{0}\left(\displaystyle \frac{1}{4!}{\eta }_{{iiii}}{x}^{4}\right){\psi }_{0}\,{dx}.\end{eqnarray} \tag{ A7 }$$

Let us begin with A first. We need to remember that when an operator involves a differentiation, it does not commute. This means we now need to first differentiate ψ₀ twice, but leave ψ₀ unaffected:

$$\begin{eqnarray}&&A={\int }_{-\infty }^{\infty }{\psi }_{0}^{* }\left(-\displaystyle \frac{1}{2}\displaystyle \frac{{\partial }^{2}{\psi }_{0}}{\partial {Q}^{2}}\right){dx}.\end{eqnarray} \tag{ A8 }$$

Now to insert the equations for the wave functions:

$$\begin{eqnarray}&&\begin{array}{l}=\ -\displaystyle \frac{1}{2}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\displaystyle \int }_{-\infty }^{\infty }{e}^{-\tfrac{\omega }{2}{\left(x-\sigma \right)}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}{e}^{-\tfrac{\omega }{2}{\left(x-\sigma \right)}^{2}}\,{dx}.\end{array}\end{eqnarray} \tag{ A9 }$$

Recalling u = x − σ, therefore, dx = du. Substitute this back in:

$$\begin{eqnarray}&&=\ -\displaystyle \frac{1}{2}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{e}^{-\tfrac{\omega }{2}{\left(u\right)}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}{e}^{-\tfrac{\omega }{2}{\left(u\right)}^{2}}\,{du}.\end{eqnarray} \tag{ A10 }$$

Doubly differentiating within the equation we get

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial x}{e}^{-\tfrac{\omega }{2}{\left(u\right)}^{2}}=-\omega {{ue}}^{-\tfrac{\omega }{2}{\left(u\right)}^{2}},\end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial x}-\omega {{ue}}^{-\tfrac{\omega }{2}{\left(u\right)}^{2}}={u}^{2}{\omega }^{2}{e}^{-\tfrac{\omega }{2}{\left(u\right)}^{2}}-\omega {e}^{-\tfrac{\omega }{2}{\left(u\right)}^{2}},\end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray}&&\therefore =-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{e}^{-\tfrac{\omega }{2}{\left(u\right)}^{2}}\left({u}^{2}{\omega }^{2}{e}^{-\tfrac{\omega }{2}{\left(u\right)}^{2}}-\omega {e}^{-\tfrac{\omega }{2}{\left(u\right)}^{2}}\right){du},\end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray}&&=\ -\displaystyle \frac{1}{2}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }\left({u}^{2}{\omega }^{2}{e}^{\omega {\left(u\right)}^{2}}-\omega {e}^{-\omega {\left(u\right)}^{2}}\right){du},\end{eqnarray} \tag{ A14 }$$

$$\begin{eqnarray}&&=\ -\displaystyle \frac{1}{2}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}\left[{\omega }^{2}\displaystyle \frac{1}{2\omega }{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2}-2\omega {\left(\displaystyle \frac{\pi }{4\omega }\right)}^{1/2}\right],\end{eqnarray} \tag{ A15 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{4}\omega .\end{eqnarray} \tag{ A16 }$$

Now onto part B, same as before we will introduce the wave functions then use substitution to solve:

$$\begin{eqnarray}&&B={\int }_{-\infty }^{\infty }{\psi }_{0}^{* }\left(\displaystyle \frac{1}{2!}{\eta }_{{ii}}{x}^{2}\right){\psi }_{0}\,{dx},\end{eqnarray} \tag{ A17 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{2!}{\eta }_{{ii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{x}^{2}{e}^{-\omega {(x-\sigma )}^{2}}\,{dx},\end{eqnarray} \tag{ A18 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{2!}{\eta }_{{ii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{(u+\sigma )}^{2}{e}^{-\omega {u}^{2}}\,{du},\end{eqnarray} \tag{ A19 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{2!}{\eta }_{{ii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }({u}^{2}+2u\sigma +{\sigma }^{2}){e}^{-\omega {u}^{2}}\,{du}.\end{eqnarray} \tag{ A20 }$$

Once again we can break this down into pieces:

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }({u}^{2}){e}^{-\omega {u}^{2}}\,{du}=\displaystyle \frac{1}{2}\left(\displaystyle \frac{\pi }{{\omega }^{3}}\right),\end{eqnarray} \tag{ A21 }$$

$$\begin{eqnarray}&&2\sigma {\int }_{-\infty }^{\infty }(u){e}^{-\omega {u}^{2}}\,{du}=0,\end{eqnarray} \tag{ A22 }$$

$$\begin{eqnarray}&&{\sigma }^{2}{\int }_{-\infty }^{\infty }{e}^{-\omega {u}^{2}}\,{du}={\sigma }^{2}\left(\displaystyle \frac{\pi }{\omega }\right).\end{eqnarray} \tag{ A23 }$$

Putting these components back together:

$$\begin{eqnarray}&&=\displaystyle \frac{1}{2!}{\eta }_{{ii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}\left(\displaystyle \frac{1}{2\omega }{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2}+{\sigma }^{2}{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2}\right),\end{eqnarray} \tag{ A24 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{1}{2!}{\eta }_{{ii}}\left[\displaystyle \frac{1}{2\omega }+{\sigma }^{2}\right].\end{eqnarray} \tag{ A25 }$$

Part C, using the same method as before:

$$\begin{eqnarray}&&C={\int }_{-\infty }^{\infty }{\psi }_{0}^{* }\left(\displaystyle \frac{1}{3!}{\eta }_{{iii}}{x}^{3}\right){\psi }_{0}\,{dx},\end{eqnarray} \tag{ A26 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{3!}{\eta }_{{iii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{x}^{3}{e}^{-\omega {(Q-\sigma )}^{2}}\,{dx},\end{eqnarray} \tag{ A27 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{3!}{\eta }_{{iii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{(u+\sigma )}^{3}{e}^{-\omega {\left(u\right)}^{2}}\,{du},\end{eqnarray} \tag{ A28 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{3!}{\eta }_{{iii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{\left({u}^{3}+3\sigma {u}^{2}+{\sigma }^{2}u+{\sigma }^{3}\right)}^{3}{e}^{-\omega {\left(u\right)}^{2}}\,du.\end{eqnarray} \tag{ A29 }$$

Breaking it down into parts again for ease of solving:

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }{u}^{3}{e}^{-\omega {\left(u\right)}^{2}}\,{du}=0,\end{eqnarray} \tag{ A30 }$$

$$\begin{eqnarray}&&3\sigma {\int }_{-\infty }^{\infty }{u}^{2}{e}^{-\omega {\left(u\right)}^{2}}\,{du}=3\sigma {\left(\displaystyle \frac{\pi }{4{\omega }^{3}}\right)}^{1/2},\end{eqnarray} \tag{ A31 }$$

$$\begin{eqnarray}&&{\sigma }^{2}{\int }_{-\infty }^{\infty }{{ue}}^{-\omega {\left(u\right)}^{2}}\,{du}=0,\end{eqnarray} \tag{ A32 }$$

$$\begin{eqnarray}&&{\sigma }^{3}{\int }_{-\infty }^{\infty }{e}^{-\omega {\left(u\right)}^{2}}\,{du}={\sigma }^{3}{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2}.\end{eqnarray} \tag{ A33 }$$

Substituting these components back into the main equation:

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{3!}{\eta }_{{iii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}\left[\displaystyle \frac{3\sigma }{2\omega }{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2}+{\sigma }^{3}{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2}\right],\end{eqnarray} \tag{ A34 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{3!}{\eta }_{{iii}}\left[\displaystyle \frac{3\sigma }{2\omega }+{\sigma }^{3}\right].\end{eqnarray} \tag{ A35 }$$

Finally, onto section D of the original equation:

$$\begin{eqnarray}&&D={\int }_{-\infty }^{\infty }{\psi }_{0}^{* }\left(\displaystyle \frac{1}{4!}{\eta }_{{iiii}}{x}^{4}\right){\psi }_{0}\,{dx},\end{eqnarray} \tag{ A36 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{4!}{\eta }_{{iiii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{x}^{4}{e}^{-\omega {(Q-\sigma )}^{2}}\,{dx},\end{eqnarray} \tag{ A37 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{4!}{\eta }_{{iiii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{(u+\sigma )}^{4}{e}^{-\omega {\left(u\right)}^{2}}\,{du},\end{eqnarray} \tag{ A38 }$$

$$\begin{eqnarray}&&\begin{array}{l}=\ \displaystyle \frac{1}{4!}{\eta }_{{iiii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}\\ \ \ \times \ {\int }_{-\infty }^{\infty }\left({u}^{4}+4\sigma {u}^{3}+6{\sigma }^{2}{u}^{2}+4{\sigma }^{3}u+{u}^{4}\right){e}^{-\omega {u}^{2}}\,{du}.\end{array}\end{eqnarray} \tag{ A39 }$$

Break down into parts:

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }{u}^{4}{e}^{-\omega {u}^{2}}{du}=\displaystyle \frac{6}{8{\omega }^{2}}{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2},\end{eqnarray} \tag{ A40 }$$

$$\begin{eqnarray}&&4\sigma {\int }_{-\infty }^{\infty }{u}^{3}{e}^{-\omega {u}^{2}}{du}=0,\end{eqnarray} \tag{ A41 }$$

$$\begin{eqnarray}&&6{\sigma }^{2}{\int }_{-\infty }^{\infty }{u}^{2}{e}^{-\omega {u}^{2}}{du}=\displaystyle \frac{6{\sigma }^{2}}{2\omega }{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2},\end{eqnarray} \tag{ A42 }$$

$$\begin{eqnarray}&&4{\sigma }^{3}{\int }_{-\infty }^{\infty }u\,{du}=0,\end{eqnarray} \tag{ A43 }$$

$$\begin{eqnarray}&&{\sigma }^{4}{\int }_{-\infty }^{\infty }{e}^{-\omega {u}^{2}}{du}={\sigma }^{4}{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2}.\end{eqnarray} \tag{ A44 }$$

Summing the pieces back together:

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{4!}{\eta }_{{iiii}}{\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}\left[\displaystyle \frac{6}{8{\omega }^{2}}{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2}+\displaystyle \frac{6{\sigma }^{2}}{2\omega }{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2}+{\sigma }^{4}{\left(\displaystyle \frac{\pi }{\omega }\right)}^{1/2}\right],\end{eqnarray} \tag{ A45 }$$

$$\begin{eqnarray}&&=\ \displaystyle \frac{1}{4!}{\eta }_{{iiii}}\left[\displaystyle \frac{3}{4{\omega }^{2}}+\displaystyle \frac{3{\sigma }^{2}}{\omega }+{\sigma }^{4}\right].\end{eqnarray} \tag{ A46 }$$

If we now take all the parts, A, B, C, and D, we can now have the full equation for the energy of the ground state, which matches the expression in Equation (6):

$$\begin{eqnarray}\begin{array}{rcl}{E}_{0} & = & \displaystyle \frac{1}{4}\omega +\displaystyle \frac{1}{2!}{\eta }_{{ii}}\left[\displaystyle \frac{1}{2\omega }+{\sigma }^{2}\right]+\displaystyle \frac{1}{3!}{\eta }_{{iii}}\left[\displaystyle \frac{3\sigma }{2\omega }+{\sigma }^{3}\right]\\ & & +\ \displaystyle \frac{1}{4!}{\eta }_{{iiii}}\left[\displaystyle \frac{3}{4{\omega }^{2}}+\displaystyle \frac{3{\sigma }^{2}}{\omega }+{\sigma }^{4}\right].\end{array}\end{eqnarray} \tag{ A47 }$$

The method described here can be applied to the first excited state to achieve the answer shown in Equation (6).

Let us first start with the expectation value for the harmonic position, u:

$$\begin{eqnarray}&&\langle u\rangle ={\int }_{-\infty }^{\infty }u| \psi {| }^{2}du.\end{eqnarray} \tag{ A48 }$$

We can use a single simple identity to calculate the expectation value, as we know the wave function can be normalized:

$$\begin{eqnarray}&&{\int }_{-t}^{t}o(u){du}=0;| \psi {| }^{2}=1,\end{eqnarray} \tag{ A49 }$$

$$\begin{eqnarray}&&\therefore \langle u\rangle =0.\end{eqnarray} \tag{ A50 }$$

Now to do the anharmonic position expectation values. Although the anharmonic wave functions are also normalized, we cannot use only the identities described previously, as the equations are more complex. Therefore, in this case we need to define the wave functions. Using Equation (2), we see that the for the ground state is

$$\begin{eqnarray}&&{\psi }_{0}={\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/4}{e}^{-\tfrac{\omega }{2}{\left(x-\sigma \right)}^{2}}.\end{eqnarray} \tag{ A51 }$$

Additionally, for the first excited vibrational state:

$$\begin{eqnarray}&&{\psi }_{1}={\left(\displaystyle \frac{4{\omega }^{3}}{\pi }\right)}^{1/4}(x-\sigma ){e}^{-\tfrac{\omega }{2}{\left(x-\sigma \right)}^{2}}.\end{eqnarray} \tag{ A52 }$$

Now to derive the expectation value for the shifted anharmonic position, 〈x〉, with respect to the anharmonic wave functions. We can now derive the expectation value of the position, x, for the anharmonic ground and first vibrational level:

$$\begin{eqnarray}&&\langle {\psi }_{0}| x| {\psi }_{0}\rangle ={\int }_{-\infty }^{\infty }| {\psi }_{0}{| }^{2}x\,{dx},\end{eqnarray} \tag{ A53 }$$

$$\begin{eqnarray}&&=\ {\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{{xe}}^{-\omega {(x-\sigma )}^{2}}\,{dx}.\end{eqnarray} \tag{ A54 }$$

Substituting x = u + σ:

$$\begin{eqnarray}&&=\ {\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }(u+\sigma ){e}^{-\omega {\left(u\right)}^{2}}\,du.\end{eqnarray} \tag{ A55 }$$

Separating into two components (A and B) for easier integration:

$$\begin{eqnarray}&&A={\int }_{-\infty }^{\infty }(u){e}^{-\omega {\left(u\right)}^{2}}\,du.\end{eqnarray} \tag{ A56 }$$

Using the identity for an odd function (where o(u) = − o( − u) for all u), and knowing the wave function can be normalized:

$$\begin{eqnarray}&&{\int }_{-t}^{t}o(u){du}=0,\end{eqnarray} \tag{ A57 }$$

$$\begin{eqnarray}&&\therefore A={\int }_{-\infty }^{\infty }(u){e}^{-\omega {\left(u\right)}^{2}}\,{du}=0.\end{eqnarray} \tag{ A58 }$$

Now for part B:

$$\begin{eqnarray}&&B=\sigma {\int }_{-\infty }^{\infty }{e}^{-\omega {\left(u\right)}^{2}}\,du.\end{eqnarray} \tag{ A59 }$$

An even function e(u) satisfies e(u) = e( − u) for all u. There, for any t,

$$\begin{eqnarray}&&{\int }_{-t}^{t}e(u){du}=2{\int }_{0}^{t}e(u)du,\end{eqnarray} \tag{ A60 }$$

$$\begin{eqnarray}&&\therefore B=\sigma {\int }_{-\infty }^{\infty }{e}^{-\omega {\left(u\right)}^{2}}\,{du}=2\sigma {\int }_{0}^{\infty }{e}^{-\omega {\left(u\right)}^{2}}\,{du}.\end{eqnarray} \tag{ A61 }$$

Then using a secondary integration identity:

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{e}^{-\omega {u}^{2}}{du}={\left(\displaystyle \frac{\pi }{4\omega }\right)}^{1/2},\end{eqnarray} \tag{ A62 }$$

$$\begin{eqnarray}&&\therefore B=2\sigma {\left(\displaystyle \frac{\pi }{4\omega }\right)}^{1/2}.\end{eqnarray} \tag{ A63 }$$

Put A and B back together into the original expression:

$$\begin{eqnarray}&&=\ {\left(\displaystyle \frac{\omega }{\pi }\right)}^{1/2}\left[0+2\sigma {\left(\displaystyle \frac{\pi }{4\omega }\right)}^{1/2}\right],\end{eqnarray} \tag{ A64 }$$

$$\begin{eqnarray}&&=\ \sigma .\end{eqnarray} \tag{ A65 }$$

Now onto the first vibrational level:

$$\begin{eqnarray}&&\langle {\psi }_{1}| x| {\psi }_{1}\rangle ={\int }_{-\infty }^{\infty }| {\psi }_{1}{| }^{2}x\,dx,\end{eqnarray} \tag{ A66 }$$

$$\begin{eqnarray}&&=\ {\left(\displaystyle \frac{4{\omega }^{3}}{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{(x-\sigma )}^{2}{{xe}}^{-\omega {(x-\sigma )}^{2}}\,dx.\end{eqnarray} \tag{ A67 }$$

Substituting x = u + σ:

$$\begin{eqnarray}&&=\ {\left(\displaystyle \frac{4{\omega }^{3}}{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }{u}^{2}(u+\sigma ){e}^{-\omega {\left(u\right)}^{2}}\,du,\end{eqnarray} \tag{ A68 }$$

$$\begin{eqnarray}&&=\ {\left(\displaystyle \frac{4{\omega }^{3}}{\pi }\right)}^{1/2}{\int }_{-\infty }^{\infty }({u}^{3}+\sigma {u}^{2}){e}^{-\omega {\left(u\right)}^{2}}\,du.\end{eqnarray} \tag{ A69 }$$

Once again separating into two components (A and B) for easier integration:

$$\begin{eqnarray}&&A={\int }_{-\infty }^{\infty }({u}^{3}){e}^{-\omega {\left(u\right)}^{2}}\,du.\end{eqnarray} \tag{ A70 }$$

Using the identity for an odd function (where o(u) = − o( − u) for all u):

$$\begin{eqnarray}&&{\int }_{-t}^{t}o(u){du}=0,\end{eqnarray} \tag{ A71 }$$

$$\begin{eqnarray}&&\therefore A={\int }_{-\infty }^{\infty }({u}^{3}){e}^{-\omega {\left(u\right)}^{2}}\,{du}=0,\end{eqnarray} \tag{ A72 }$$

$$\begin{eqnarray}&&B=\sigma {\int }_{-\infty }^{\infty }({u}^{2}){e}^{-\alpha {\left(u\right)}^{2}}\,{du}.\end{eqnarray} \tag{ A73 }$$

Using the integration identity:

$$\begin{eqnarray}&&{\int }_{-\infty }^{\infty }{u}^{2}{e}^{-\omega {u}^{2}}{du}=\displaystyle \frac{1}{2}{\left(\displaystyle \frac{\pi }{{\omega }^{3}}\right)}^{1/2},\end{eqnarray} \tag{ A74 }$$

$$\begin{eqnarray}&&\therefore B=\displaystyle \frac{\sigma }{2}{\left(\displaystyle \frac{\pi }{{\omega }^{3}}\right)}^{1/2}.\end{eqnarray} \tag{ A75 }$$

Putting A and B back together into the original expression:

$$\begin{eqnarray}&&=\ {\left(\displaystyle \frac{4{\omega }^{3}}{\pi }\right)}^{1/2}\left[0+\displaystyle \frac{\sigma }{2}{\left(\displaystyle \frac{\pi }{{\omega }^{3}}\right)}^{1/2}\right],\end{eqnarray} \tag{ A76 }$$

$$\begin{eqnarray}&&=\ \sigma .\end{eqnarray} \tag{ A77 }$$

Clearly the outcome is independent of the vibrational level.

This next section will detail the working out required to produce Equations (26), (27), (28), (29) and (30). The general analytical expression for a spectrum of a diatomic molecule is the following:

$$\begin{eqnarray}&&{\varepsilon }_{(J,v)}={\varepsilon }_{J}+{\varepsilon }_{v},\end{eqnarray} \tag{ A78 }$$

$$\begin{eqnarray}&&=\ {\bar{\omega }}_{v}(v+1/2)+{B}_{v}J(J+1)-{D}_{e}{J}^{2}{\left(J+1\right)}^{2},\end{eqnarray} \tag{ A79 }$$

where

$$\begin{eqnarray}&&{\bar{\omega }}_{v}=\ {\bar{\omega }}_{e}-{\bar{\omega }}_{e}{x}_{e}(v+1/2),\end{eqnarray} \tag{ A80 }$$

$$\begin{eqnarray}\begin{array}{rcl}\therefore {\varepsilon }_{(J,v)} & = & {\bar{\omega }}_{e}(v+1/2)-{\bar{\omega }}_{e}{x}_{e}{\left(v+1/2\right)}^{2}\\ & & +\ {B}_{v}J(J+1)-{D}_{e}{J}^{2}{\left(J+1\right)}^{2}+....\end{array}\end{eqnarray} \tag{ A81 }$$

The frequency ν_0→1 is usually called the band origin or band center, which is represented by the following equation:

$$\begin{eqnarray}&&{\nu }_{0\to 1}={\bar{\omega }}_{e}-2{\bar{\omega }}_{e}{x}_{e}.\end{eqnarray} \tag{ A82 }$$

D_e represents the centrifugal distortion constant. For current theory we can ignore the small centrifugal distortions from D_e , etc., and the equation can be rewritten simply as

$$\begin{eqnarray}&&{\varepsilon }_{(J,v)}={\bar{\omega }}_{e}(v+1/2)-{\bar{\omega }}_{e}{x}_{e}{\left(v+1/2\right)}^{2}+{B}_{v}J(J+1).\end{eqnarray} \tag{ A83 }$$

For the next part we will be restricting the discussion to solely fundamental vibrational transitions, v = 0 → v = 1. We also take the respective B values as B₀ and B₁ with B₀ > B₁. This transition can generally be described as (in cm⁻¹):

$$\begin{eqnarray}&&{\rm{\Delta }}\varepsilon ={\varepsilon }_{J^{\prime} ,v=1}-{\varepsilon }_{J^{\prime\prime} ,v=0},\end{eqnarray} \tag{ A84 }$$

$$\begin{eqnarray}&&=\ {\nu }_{0\to 1}+{B}_{1}J^{\prime} (J^{\prime} +1)-{B}_{0}J^{\prime\prime} (J^{\prime\prime} +1).\end{eqnarray} \tag{ A85 }$$

The transitions for the P branch (ΔJ = − 1, $J^{\prime\prime} =J^{\prime} +1$) is given by

$$\begin{eqnarray}&&{\rm{\Delta }}\varepsilon ={\bar{\nu }}_{P},\end{eqnarray} \tag{ A86 }$$

$$\begin{eqnarray}&&=\ {\nu }_{0\to 1}-({B}_{1}+{B}_{0})(J^{\prime} +1)+({B}_{1}-{B}_{0}){\left(J^{\prime} +1\right)}^{2},\end{eqnarray} \tag{ A87 }$$

where $J^{\prime} =0,1,2,3,..$. The transitions for the Q branch (ΔJ = 0, $J^{\prime\prime} =J^{\prime}$) is given by

$$\begin{eqnarray}&&{\rm{\Delta }}\varepsilon ={\bar{\nu }}_{Q},\end{eqnarray} \tag{ A88 }$$

$$\begin{eqnarray}&&=\ {\nu }_{0\to 1}-J^{\prime} (J^{\prime} +1)({B}_{1}-{B}_{0}),\end{eqnarray} \tag{ A89 }$$

where $J^{\prime} =0,1,2,3,..$.

Therefore for the R branch (ΔJ = + 1, $J^{\prime} =J^{\prime\prime} +1$) the transition is given as

$$\begin{eqnarray}&&{\rm{\Delta }}\varepsilon ={\bar{\nu }}_{R},\end{eqnarray} \tag{ A90 }$$

$$\begin{eqnarray}&&=\ {\nu }_{0\to 1}+({B}_{1}+{B}_{0})(J^{\prime\prime} +1)+({B}_{1}-{B}_{0}){\left(J^{\prime\prime} +1\right)}^{2},\end{eqnarray} \tag{ A91 }$$

where J'' = 0, 1, 2, 3, ....

Taking the general rovibrational equation:

$$\begin{eqnarray}&&\varepsilon (J,v)={\bar{\omega }}_{e}(v+1/2)-\bar{\omega }{x}_{e}{\left(v+1/2\right)}^{2}+{B}_{v}J(J+1).\end{eqnarray} \tag{ A92 }$$

Let us now apply this for the ground state with quantum number J'' and the upper vibrational state with quantum number J' as per standard notation:

$$\begin{eqnarray}&&{\varepsilon }_{(J^{\prime\prime} ,v=0)}=\displaystyle \frac{1}{2}{\bar{\omega }}_{e}-\displaystyle \frac{1}{4}\bar{\omega }{x}_{e}+{B}_{0}J^{\prime\prime} (J^{\prime\prime} +1),\end{eqnarray} \tag{ A93 }$$

$$\begin{eqnarray}&&{\varepsilon }_{(J^{\prime} ,v=1)}=\displaystyle \frac{3}{2}{\bar{\omega }}_{e}-\displaystyle \frac{9}{4}\bar{\omega }{x}_{e}+{B}_{1}J^{\prime} (J^{\prime} +1),\end{eqnarray} \tag{ A94 }$$

$$\begin{eqnarray}&&\therefore {\rm{\Delta }}\varepsilon ={\nu }_{0\to 1}+{B}_{1}J^{\prime} (J^{\prime} +1)-{B}_{0}J^{\prime\prime} (J^{\prime\prime} +1).\end{eqnarray} \tag{ A95 }$$

Now consider the S branch where ${\rm{\Delta }}J=+2,J^{\prime} =J^{\prime\prime} +2$. By plugging in these formulas:

$$\begin{eqnarray}&&{\rm{\Delta }}\varepsilon ={\nu }_{S},\end{eqnarray} \tag{ A96 }$$

$$\begin{eqnarray}&&=\ {\nu }_{0\to 1}+{B}_{1}(J^{\prime\prime} +2)(J^{\prime\prime} +3)-{B}_{0}J^{\prime\prime} (J^{\prime\prime} +2),\end{eqnarray} \tag{ A97 }$$

$$\begin{eqnarray}&&=\ {\nu }_{0\to 1}+{B}_{1}({\left(J^{\prime\prime} \right)}^{2}+5J^{\prime\prime} +6)-{B}_{0}({\left(J^{\prime\prime} \right)}^{2}+J^{\prime\prime} ).\end{eqnarray} \tag{ A98 }$$

In the cases where B₀ = B₁ = B, the expression simplifies to

$$\begin{eqnarray}&&{\nu }_{S}={\nu }_{0\to 1}+B(4J^{\prime\prime} +6).\end{eqnarray} \tag{ A99 }$$

Finally, onto the O branch where ${\rm{\Delta }}J=-2,J^{\prime\prime} =J^{\prime} +2$. By plugging in these formulas:

$$\begin{eqnarray}&&{\rm{\Delta }}\varepsilon ={\nu }_{O},\end{eqnarray} \tag{ A100 }$$

$$\begin{eqnarray}&&=\ {\nu }_{0\to 1}+{B}_{1}J^{\prime} (J^{\prime} +1)-{B}_{0}(J^{\prime} +2)(J^{\prime} +3),\end{eqnarray} \tag{ A101 }$$

$$\begin{eqnarray}&&=\ {\nu }_{0\to 1}+{B}_{1}{({(J^{\prime} )}^{2}+J^{\prime} )-{B}_{0}((J^{\prime} )}^{2}+5J^{\prime} +6).\end{eqnarray} \tag{ A102 }$$

In the cases where B₀ = B₁ = B, the expression simplifies to

$$\begin{eqnarray}&&{\nu }_{O}=\ {\nu }_{0\to 1}-B(4J^{\prime\prime} +6).\end{eqnarray} \tag{ A103 }$$