An ab initio study of the rovibronic spectrum of sulphur monoxide (SO): diabatic vs. adiabatic representation

Abstract

We present an ab initio study of the rovibronic spectra of sulphur monoxide (³²S¹⁶O) using internally contracted multireference configuration interaction (ic-MRCI) method and aug-cc-pV5Z basis sets. It covers 13 electronic states X³Σ⁻, a¹Δ, b¹Σ⁺, c¹Σ⁻, A′′³Σ⁺, A′³Δ, A³Π, B³Σ⁻, C³Π, d¹Π, e¹Π, C′³Π, and (3)¹Π ranging up to 66 800 cm⁻¹. The ab initio spectroscopic model includes 13 potential energy curves, 23 dipole and transition dipole moment curves, 23 spin–orbit curves, and 14 electronic angular momentum curves. A diabatic representation is built by removing the avoided crossings between the spatially degenerate pairs C³Π–C′³Π and e¹Π–(3)¹Π through a property-based diabatisation method. We also present non-adiabatic couplings and diabatic couplings for these avoided crossing systems. All phases for our coupling curves are defined, and consistent, providing the first fully reproducible spectroscopic model of SO covering the wavelength range longer than 147 nm. Finally, an ab initio rovibronic spectrum of SO is computed.

---

1 Introduction

Sulphur monoxide (SO) has many experimental and observational applications, making it a molecule of continuing interest in spectroscopic and kinetic studies. SO is an important intermediate in the oxidation of sulphur compounds within combustion reactions,¹ making SO studies applicable to environmental chemistry since sulphur oxides lead to the acidification/pollution of Earth's atmosphere.^2–5 Burkholder et al.⁶ show the importance of SO within chemical reactions of Earth's atmosphere due to its high reactivity and chemical involvement with N₂ and O₂. Additionally, UV lasing in SO was demonstrated through optically pumping transitions within the B³Σ⁻–X³Σ⁻ electronic band⁷ and then later for UV energy storage lasers.⁸ SO also has astrophysical importance, its Zeeman splitting has been used to probe the magnetic fields of dense star-forming regions, such as the Orion molecular cloud, for field strengths ≥10⁻² G^9,10 and its presence within warm chemistry's mean it is an excellent shock tracer.^11,12

The frontier orbitals of SO resemble that of carbon monoxide where the two π* valence electrons mean SO has a X³Σ⁻ ground state similar to O₂ and S₂. Being isoelectronic with O₂, SO has two low-lying metastable a¹Δ and b¹Σ⁺ states which are relatively short lived due to large spin–orbit coupling. SO has been subject to pure rotational,^13–15 electronic,^16–20 and ro-vibrational^21–23 experimental spectroscopic studies. More recently, Bernath et al.^24,25 produced empirical line lists for the b¹Σ⁺–X³Σ⁻ and a¹Δ–X³Σ⁻ rovibronic bands and for the X³Σ⁻ and a¹Δ rovibrational bands of SO. The literature also contains many theoretical studies of SO, the most comprehensive, and one we compare to often, being Yu and Bian²⁶ who give spectroscopic constants on all electronic states we consider here (except (3)¹Π) computed through internally contracted multireference configuration interaction (ic-MRCI) methods using aug-cc-pV5Z basis sets. Another important theoretical work is by Sarka and Nanbu²⁷ who study the UV region of SO non-adiabatically where they compute PECs, DMCs, and non-adiabatic couplings (NACs) for the X³Σ⁻, A³Π, B³Σ⁻, C³Π, C′³Π, (3)³Σ⁻, (4)³Π, and (5)³Π states at a MRCI-F12+Q level of theory using aug-cc-pV(5+d)Z basis sets. Sarka and Nanbu²⁷ are also the first to compute cross-sections for SO longer than 190 nm in the UV.

Despite the strong experimental, theoretical, and observational baseline for SO, to this date the existing spectroscopic line list data for SO are limited in coverage. HITRAN²⁸ contains only transitions between electronic states X³Σ⁻, a¹Δ and b¹Σ⁺ considering relatively low vibrational excitations. SO data and spectroscopic databases NIST²⁹ and CDMS³⁰ contain data covering the microwave region only. In this study we aim to provide a strong theoretical base from which an accurate spectroscopic description of SO with large coverage can be made through the refinement of our ab initio model to empirically determined energies in a future work. To achieve this we calculate ab initio potential energy curves (PECs), spin–orbit curves (SOCs), electronic angular momentum curves (EAMCs), and electric (transition) dipole moment curves ((T)DMCs) for the 13 lowest electronic states of SO (X³Σ⁻, a¹Δ, b¹Σ⁺, c¹Σ⁻, A′³Δ, A′′³Σ⁺, A³Π, B³Σ⁻, C³Π, d¹Π, e¹Π, C′³Π, (3)¹Π) at an MRCI level of theory using aug-cc-pV5Z basis sets. The relative phases of the the non-diagonal couplings and transition dipole moments provided are fully self-consistent,³¹ which is crucial for reproducible spectroscopic studies.

Our ab initio curves of SO are adiabatic as computed under the Born–Oppenheimer approximation³² and so the spatially degenerate states e¹Π, (3)¹Π and C³Π, C′³Π of SO exhibit avoided crossings due to their shared symmetries, where the non-adiabatic effects play important role. The associated non-adiabatic couplings (NACs) arise from the nuclear motion kinetic energy operator acting on electronic wavefunctions and contain first- and second-order differential operators with respect to the nuclear coordinates. Inclusion of NACs within numerical dynamics calculations are computationally costly when using analytical forms for the PECs and couplings because of both the cusp-like nature of PECs and the NACs having singularities at the region of spatial degeneracy.^33–36 This makes the fitting of analytical forms to the PECs and couplings almost impossible. Here we explore a diabatisation procedure to transform the adiabatic basis to the diabatic basis using NACs reconstructed from PECs as opposed to obtaining them ab initio,²⁷ where the adiabatic first-order differential couplings are negligible or vanish, at the cost of inducing off-diagonal potential energy couplings,³⁷ or diabatic couplings (DCs).^33,34 In this representation, the electronic wavefunctions are smooth which lessens both numerical problems encountered in calculations within the adiabatic representation and the computational cost in computing NACs. Diabatisation helps recover non-Born–Oppenheimer dynamics^34,38 and is important for processes such as collisions between open-shell molecules with spatially degenerate electronic states.^35,39–42 The importance of the non-adiabatic effects on the spectral properties of SO is analysed.

2 Computational details

Internally-contracted multi-reference configuration interaction (icMRCI) ab initio calculations for the 13 lowest states of SO correlating with S(³P) + O(³P), S(¹D) + O(³P) and S(¹D) + O(¹D) were performed using Molpro⁴³ with aug-cc-pV5Z basis sets,^44,45 using molecular orbitals obtained from state-averaged complete active space self-consistent field (CASSCF) calculations. Under C_2V point group symmetry all ab initio calculations started with 14 (8σ, 3π_x, 3π_y) orbitals which included 6 closed (4σ, 1π_x, 1π_y) orbitals. The active space occupying 12 active electrons consisted of 8 (5σ–8σ, 2π, 3π) valence orbitals. The PECs, including the 8 bound states X³Σ⁻, a¹Δ, b¹Σ⁺, c¹Σ⁻, A′′³Σ⁺, A′³Δ, A³Π, and B³Σ⁻ are shown in Fig. 1, as well as the adiabatic e¹Π–(3)¹Π and C³Π–C′³Π systems. The EAMC, SOCs (both diagonal and non-diagonal), DMCs (diagonal and transition) computed and used later for cross-section calculations are shown in Fig. 2–6 in the original adiabatic representation as computed by Molpro 2015⁴³ (left) and the diabatic representation (right). Further discussion of the diabatisation is given in the next section.

Fig. 1 Plots of our 13 ab initio PECs covering the first 7 triplet and 6 singlet electronic states up to 80 000 cm⁻¹ over inter-nuclear separations 1.0–4.5 Å. The adiabatic and diabatic representations of the PECs are shown in the left and right hand panels respectively.

Fig. 2 Ab initio spin–orbit curves between states of the same Λ (projection of the angular momentum) as a function of bond length. In general, the spin–orbit coupling is strong also between states of different multiplicity. These curves have been multiplied by −i for visual purposes.

Fig. 3 Ab initio spin–orbit MOLRPO matrix elements in the adiabatic (left) and diabatic (right) representations between states of different values of Λ as a function of bond length. These curves have been multiplied by −i for visual purposes.

Fig. 4 Ab initio (transition) dipole moment matrix elements (Debye) between states of the same symmetry (Λ and multiplicity) in the adiabatic (left) and diabatic (right) representations as a function of the bond length (Å).

Fig. 5 Ab initio transition dipole moment curves (Debye) between states of different symmetry in the adiabatic (left) and diabatic (right) representations as a function of the bond length (Å).

Fig. 6 Ab initio electronic angular momentum curves in the adiabatic (left) and diabatic (right) representation plotted over bond length. L_x means the x-component of the electronic angular momentum. These curves have been multiplied by (−i) for visual purposes.

3 Diabatisation

3.1 Non-adiabatic couplings (NACs)

The commonly employed Born–Oppenheimer approximation assumes the nuclear dynamics evolve adiabatically on a single, uncoupled electronic potential energy surface. This results in the appearance of avoided crossings between electronic states that are degenerate in energy at some geometry and large gradients in couplings connecting these states. The Born–Huang expansion⁴⁶ allows one to recover the non-adiabatic behaviour in the region of the avoided crossing by including non-adiabatic couplings (NACs) which introduce off-diagonal kinetic energy terms in the molecular Hamiltonian. An alternative approach is to transform the electronic Hamiltonian to the diabatic basis, where electronic states are coupled via off-diagonal potential energy terms, termed diabatic couplings (DCs).^33,34,37 The transformation from the adiabatic to the diabatic basis is described by a unitary matrix U, which is parametrically dependent on the NAC term,where the mixing angle β(R) is obtained by integrating the functional form of the non-adiabatic derivative coupling ^35,47–49with |ψ₁〉 and |ψ₂〉 as the lower and upper energy electronic wavefunctions in the adiabatic representation.

Writing the two-state electronic Hamiltonian in terms of the adiabatic potential energy curves V^a₁(R) and V^a₂(R),

we can obtain the diabatic Hamiltonian by applying the unitary transformation U(β(R)),where the superscripts ‘d’ and ‘a’ refer to the diabatic and adiabatic bases, respectively, and the off-diagonal elements V^d₁₂(R) are the DCs. The reverse transformation is obtained by diagonalising the diabatic representation of the electronic Hamiltonian.

The NAC can be computed via quantum-chemistry methods from the electronic wavefunctions, as done by Sarka and Nanbu²⁷ for SO. It is symmetrical with a cusp at the crossing point R_c. Alternatively, the NAC curves are modelled using, e.g. a Lorentzian,⁵⁰ as given by

where α is the inverse half-width-at-half-maximum (HWHM), or a Laplacianwhere γ is a damping constant related to the HWHM, superscripts ‘Lo’ and ‘La’ mean Lorentzian and Laplacian respectively, and the normalisation is applied. Fig. 7 illustrates the C³Π–C′³Π NAC modelled in this work using a Lorentzian and Laplacian function. The mixing angle β(R), determined through eqn (2), ranges from 0 to π/2 going through π/4 at the crossing point R = R_c can also be seen in Fig. 7.

Fig. 7 Comparison of example NACs (ϕ_ij) and corresponding diabatic mixing angles (β) between the Lorentzian (‘Lo’), Laplacian (‘La’), and geometrically averaged (‘ga’) cases as described in the text. These curves are computed for the C³Π and C′³Π non-adiabatic interaction (see Section 3.2.1 and Fig. 1).

The Lorentzian was shown to provide a good description of the ab initio NACs around the crossing point^48,49,51,52 (see Fig. 7) but diverges at large distances R from R_c causing improper shaped diabatic PECs by decaying too slowly.^48,49,51,53 It has been discussed that some damping functions can be introduced to correct the Lorentzian's slow decay using properties such as dipole moments, but determination of their fitting parameters is both difficult and requires extra calculations.^50,51 Laplacians underestimate NACs far in the wings and overestimate them near the crossing point.⁴⁸ One can show that NACs have an overlap dependence on the internuclear separation,^51,54R, which is given properly by a Laplacian.⁴⁸ The undesirable features of these NAC models can be mitigated by their combination^48,49,51,53 to which we base our diabatisation procedure on. Our method of augmenting the Lorentzian with a Laplacian is discussed in Section 3.2.

3.2 Diabatisation

Here we explore the so-called ‘property based diabatisation’ method³⁵ and construct diabatic potentials using the condition of having no avoiding crossing, which we define as to minimise their second derivatives in the vicinity of the crossing point R_c:hence creating the smoothest PECs V^d₁(R) and V^d₂(R).

In order to mitigate the undesirable properties of the Lorentzian and Laplacian functional forms, we follow the approach by An and Baeck⁴⁸ and represent the mixing angle β by the following combination of the mixing angles determined from the Lorentzian and Laplacian NACs in eqn (2), β^Lo(R) and β^La(R)

where the ‘ga’ superscript refers to the geometrically averaged diabatic mixing angle (see Fig. 7). Eqn (8) must not be taken as the geometric average of β^Lo and β^La, but rather originates from the simple geometric average of V^Lo₁₂ and V^La₁₂ (see An and Baeck⁴⁸ and Appendix). An and Baeck⁴⁸ also showed that an optimal relation between the parameters α and γ exists which given by

α × γ = 1.397 (9)

providing maximal overlap between the Lorentzian and Laplacian functions over the bond length.

Where our method diverges from that of An and Baeck⁴⁸ is in the determination of the crossing point R_c and the Lorentzian parameter α. An and Baeck⁴⁸ obtained R_c and α through fitting a Lorentzian to a NAC computed with Molpro. Instead, we determine a set of parameters {R_c,α} by fulfilling the condition given in eqn (7), to which the Laplacian parameter γ is instantly obtained through eqn (9). Using the theory developed in Section 3.1 and eqn (8) the diabatising transformation U^ga corresponding to the ‘geometrically averaged’ NAC is found. With this the diabatic potential energies and DC elements can be obtained through the simple matrix transformation in eqn (4). The diabatic PECs for SO can be seen in Fig. 1 and a closeup of the avoided crossings between the e¹Π–(3)¹Π and C³Π–C′³Π states of SO superimposed with their DCs, V^ga_ij, and NACs, ϕ^ga_ij, are illustrated in Fig. 8. Fig. 8 shows that the pair of singlet states are coupled more strongly than the triplet states, and also reveals the DCs to be slightly asymmetric. This is to be expected since the DCs depend on the difference V^a₂–V^a₁ which can be asymmetrical in nature. We see an effect of this especially in the DC between the triplet states where the adiabatic PEC turning points are offset to each other by ∼0.01 Å.

Fig. 8 Illustration of the avoided crossings between e¹Π–(3)¹Π and (right panels) and C³Π–C′³Π states of SO (left panels) are shown, where adiabatic PECs are presented as dashed lines and diabatic ones in solid lines. Superimposed are the DCs (V^ga₁₂, bottom panels), and NACs (ϕ^ga_ij, top panels).

3.2.1 Potential energy curves. Fig. 1 presents ab initio PECs of the 13 lowest energy electronic states of SO. The C³Π exhibits an avoided crossing at R ∼ 2.05 Å due to a non-adiabatic coupling with the dissociative C′³Π, which lends a dissociative character to the C³Π in the long-range region. Similarly, the e¹Π state exhibits an avoided crossing at 1.95 Å, due to the singlet state (3)¹Π.²⁶ These non-adiabatic interactions produce large gradients in coupling curves connecting these states within the region of the avoided crossing, as shown for EAMCs, SOC, and DMCs in Fig. 2–6. The equilibrium geometry of the C³Π state also lies very close in energy to those of the d¹Π and B³Σ⁻ states, and so we can expect perturbations in the energy levels around their minima which was reported and confirmed experimentally by Liu et al.⁵⁵ It is worth noting that the A′³Δ and A′′³Σ⁺ states lie very close across the range of nuclear geometries included in these calculations, but do not cross.

Table 1 provides values for the optimized NAC parameters α,γ,R_c used to diabatise the energy degenerate pairs, which are visualised in Fig. 8. We see that the effect of diabatisation is to smooth the PECs, as enforced by eqn (7) with no avoided-crossing. The non-Born–Oppenheimer dynamics, initially manifested in the nuclear kintetic energy, has been rotated into the potential and coupling terms connecting the non-adiabatically interacting states. We now produce PECs that can be easily modelled by analytical forms in the diabatic representation, allowing us to tune them to experimental data in a future study where we aim to produce an empirically accurate line list.

Table 1 Optimised parameters α (inverse Lorentzian HWHM), γ (Laplacian damping parameter), and R_c (avoided crossing/centroid position) for the Lorentzian and Laplacian NACs used to diabatise the e¹Π and C³Π state PECs (see Fig. 1)

State	α	γ	R c (Å)
e^1Π	52.422	0.027	1.949
C^3Π	39.859	0.035	2.047

---

Table 2 compares the equilibrium potential energies, T_e (cm⁻¹), equilibrium bond lengths, R_e (Å), and dissociation energies, D_e (eV), of the 11 lowest singlet and triplet states of SO determined directly from the ab initio adiabatic PECs presented in Fig. 1 to both calculations^{26,27,56–58} and experiment.^{24,57,59–61} Our bond lengths show good agreement to both theoretical and empirical values, with better agreement to experiment than previous calculations for the b¹Σ⁺, A³Π, and d¹Π states. Worse agreements are seen for our T_e values, the most accurate being T_e(A³Π), T_e(C³Π), and T_e(e¹Π) within 152 cm⁻¹, 210 cm⁻¹, and 338 cm⁻¹ of experiment respectively. Lastly, we see worse agreements between our ab initio dissociation energies to both experiment and calculations, to which we underestimate. It is apparent that our dissociation asymptotes are the major source of inaccuracy in our ab initio SO model and probably arise because we do not include Sulfur specific diffuse functions in our basis set during ab initio calculations. A future goal will be to mitigate the undesirable PEC features by refining our ab initio model to experimental transition data, to which we will address when producing the final SO line list. We note that the reported Bernath et al.²⁴R_e values in Table 2 were derived from the B_ν=0 rotational constant, and show close agreements to within 0.006 cm⁻¹ Å and 0.002 cm⁻¹ Å to our bond lengths for the a¹Δ and b¹Σ⁺ respectively. Our fundamental vibrational energy of the X³Σ⁻ state is found approximately 30 cm⁻¹ too high from the experiment, which is to be expected with MRCI calculations, and provides insight to the accuracy of the other computed states and couplings, which will also require empirical tuning.

Table 2 Comparison of ab initio values of the equilibrium potential energies T_e (cm⁻¹), the dissociation energy D_e (eV) to the adiabatically correlated asymptotes, and equilibrium bond length R_e (Å) from this work to the values from the literature. Parameters next to bold state symbols correspond to the ab initio PECs calculated in this study. Calculations^26,56–58 use the MRCI/cc-aug-pV5Z level of theory, calculation by Sarka and Nanbu²⁷ use the MRCI-F12+Q/cc-aug-pV(5+d)Z level of theory, and experiments are photoion-photoelectron coincidence,⁶⁰ multiphoton ionization⁵⁷ and Ar + SO₂ afterglow,⁶¹ spectroscopes as well as from the recent analysis by Bernath et al.²⁴

State	T e (cm^−1)	D e (eV)	R e (Å)	State	T e (cm^−1)	D e (eV)	R e (Å)
X^3Σ^−	0	5.1253	1.4821	A′^3Δ	29 097.8878	1.503	1.7571
Calc.^26	0	5.418	1.4865	Calc.^26	29 828	1.72	1.7649
Expt.^59		5.429	1.481	B^3Σ^−	43 255.0097	0.9212	1.8121
Calc.^27	0	5.475	1.4925	Calc.^27	41 706.5886	1.4022	1.7868
Expt.^60			1.481	Calc.^26	41 314	1.387	1.782
Calc.^56			1.481	Expt.^60	41 629	1.410^62	1.775
Calc.^58			1.493	Calc.^58	41 206		1.794
a^1Δ	5479.8013	4.4486	1.4979	d^1Π	45 309.0766	0.0587	1.545
Calc.^26	5936	4.682	1.4945	Calc.^26	44 166	0.189	1.5475
Expt.^59	5873	4.647	1.4919	Expt.^57	43 902	0.195	1.5303
Calc.^58	5883		1.502	Calc.^58	44 975	0.059	1.723
Expt.^24			1.4920	Calc.^57	44 471	0.14	1.553
b^1Σ^+	9774.1938	3.9154	1.5057	A′′^3Σ^+	29 731.2077	1.4417	1.765
Calc.^26	10 548	4.112	1.5062	Calc.^26	30 495	1.637	1.7701
Expt.^59	10 510	4.137	1.5001	Expt.^61	30 692
Calc.^58	10 576		1.514	Calc.^58	30 025		1.776
Expt.^24			1.5035	C^3Π	44 719.2593	0.5489	1.6786
A^3P	38 607.6737	0.4246	1.6079	Calc.^26	44 033	0.609	1.6692
Calc.^27	38 879.2948	0.6441	1.5946	Calc.^27	44 909.0901	0.6027	1.6727
Calc.^26	38 334	0.665	1.6196	Expt.^57	44 381		1.654
Expt.^59	38 455	0.662	1.609	Calc.^58	44 038		1.681
Calc.^56	38 880		1.613	e^1Π	51 347.9346	0.4524	1.6864
Calc.^58	38 931		1.719	Calc.^26	51 224	0.45	1.6826
c^1Σ^−	27 274.9752	1.7679	1.7571	Expt.^57	51 558		1.6774
Calc.^26	28 544	1.879	1.7617	Calc.^58	51 258		1.685

---

The intersections between states of different symmetries obtained in our calculations can be seen in Fig. 9. The d¹Π and C³Π states cross at 1.59 and 1.82 Å, d¹Π and B³Σ⁻ at 1.64 and 2.20 Å, both agreeing with Yu and Bian²⁶ intersection locations of 1.62 and 1.80 Å and 1.60 and 2.14 Å, respectively. The intersection of the e¹Π with the B³Σ⁻ state occurs at 2.4 Å in our calculations, somewhat larger than the value of 2.3 Å reported by Yu and Bian.²⁶ Since the e¹Π and d¹Π states become repulsive at 1.92 Å and 1.9 Å respectively, crossings beyond these geometries provide potential predissociation pathways for the C³Π and B³Σ⁻ states. Yu and Bian²⁶ also show that the C′³Π state crosses the B³Σ⁻ state at 2.25 Å, to which they state predissociation of B³Σ⁻ through this C′³Π state is possible. Sarka and Nanbu²⁷ also give intersections R(C,B) = 1.57,2.21 Å as opposed to our values of R(C,B) = 1.67, 2.18 Å.

Fig. 9 Our ab initio PECs in the region 37 000–57 000 cm⁻¹ showcasing the various state crossings.

Lastly, we report further crossings between the c¹Σ⁻, A′³Δ, and A′′³Σ⁺ states and the A³Π and d¹Π states of SO at R(c,A) = 1.46 Å, R(A′,A) = 1.48 Å, R(A′′,A) = 1.50 Å, R(c,d) = 1.40 Å, R(A′,d) = 1.42 Å, and R(A′′,d) = 1.43 Å. These crossings agree with those reported by Yu and Bian²⁶ between c¹Σ⁻, A′³Δ, and A′′³Σ⁺ with A³Π in the region 1.47–1.51 Å and crossings between c¹Σ⁻, A′³Δ, A′′³Σ⁺ with d¹Π in the region 1.42–1.45 Å.

3.2.2 Spin–orbit and electronic angular momentum curves. Fig. 6 shows the EAMCs of SO both in the adiabatic and diabatic representations, where their relative phases are carefully maintained according with their ab initio values.³¹ Without this, using couplings and any other non-diagonal properties in rovibronic calculations become meaningless. The phase of the EAMCs often changes after the crossing point, lending different long-range total angular momenta of the S + O atomic system important in dissociation mechanisms.

Fig. 2 and 3 plot the z (SO_z) and x (SO_x) components of the spin–orbit curves of SO over nuclear geometries where the former couple states of same values of Λ (projection of the electronic angular momentum) and the latter couple states of different Λ. We see again that the effect of diabatisation is to smooth out the curves over R, where avoided crossings in the adiabatic picture create steep gradients around the avoided crossing. An example of the diabatisation process can be seen for the 〈e¹Π|SO_x|X³Σ⁻〉 and 〈(3)¹Π|SO_x|X³Σ⁻〉 pair in Fig. 10. Spin–orbit matrix elements at the internuclear separation of the PEC crossings are important in determining the possible spin–orbit induced predissociation mechanisms that occur between states of different spin multiplicity²⁶ (see Discussion below). Referring to the spin–orbit couplings 〈e¹Π|SO_x|B³Σ⁻〉 and 〈d¹Π|SO_x|B³Σ⁻〉 in Fig. 3 with the magnitudes of approximately 90 and 20–30 cm⁻¹, respectively, the predissociation of the B³Σ⁻ state through d¹Π is likely to be very weak, but will be stronger through the e¹Π state. The construction of diabatic SOCs will hence influence the efficiency of pre-dissociation pathways between states of different symmetry (Table 3).

Fig. 10 The diabatisation process for the curves 〈e¹Π|SO_x|X³Σ⁻〉 and 〈(3)¹Π|SO_x|X³Σ⁻〉, which become smooth as a result of the diabatisation. These curves have been multiplied by −i for visual purposes.

Table 3 Σ values (total spin angular momentum projection onto the internuclear axis) for the bra and ket electronic states of the SOCs presented in Fig. 2 and 3

SO coupling	bra Σ	ket Σ	SO coupling	bra Σ	ket Σ	SO coupling	bra Σ	ket Σ
〈A^3Π|SOx|X^3Σ^−〉	0	1	〈e^1Π|SOx|X^3Σ^−〉	0	1	〈b^1Σ^+|SOz|B^3Σ^−〉	0	0
〈A^3Π|SOx|A′^3Δ〉	0	1	〈e^1Π|SOx|B^3Σ^−〉	0	1	〈A′^3Δ|SOz|A′^3Δ〉	1	1
〈C^3Π|SOx|A′^3Δ〉	0	1	〈C^3Π|SOx|A′′^3Σ^+〉	0	1	〈A′′^3Σ^+|SOz|X^3Σ^−〉	1	1
〈A^3Π|SOx|B^3Σ^−〉	0	1	〈e^1Π|SOx|A′^3Δ〉	0	1	〈C^3Π|SOz|A^3Π〉	1	1
〈C^3Π|SOx|X^3Σ^−〉	0	1	〈d^1Π|SOx|X^3Σ^−〉	0	1	〈A^3Π|SOz|A^3Π〉	1	1
〈C^3Π|SOx|B^3Σ^−〉	0	1	〈d^1Π|SOx|B^3Σ^−〉	0	1	〈A′′^3Σ^+|SOz|B^3Σ^−〉	1	1
〈A^3Π|SOx|A′′^3Σ^+〉	0	1	〈b^1Σ^+|SOz|X^3Σ^−〉	0	0	〈C^3Π|SOz|C^3Π〉	1	1
〈d^1Π|SOx|A′^3Δ〉	0	1	〈a^1Δ|SOz|A′^3Δ〉	0	0

---

3.2.3 Dipole moment curves. Fig. 4 and 5 plot the z- and y-components of the dipole moments coupling states of same and different symmetry (Λ plus multiplicity) respectively. The corresponding μ_x components can be always reconstructed from μ_y using their symmetry properties. We see again that the effect of diabatisation is to smooth out the curves over R, where now the DMCs tend to zero in the long range limit with no steep gradients caused by avoided crossings.

The vibronic intensities are directly affected by the derivatives of the dipole moment with respect to the internucelar separation, R. Since adiabatic curves are prone to strong, steep-gradient variations around avoided crossings, even small inaccuracies in ab initio calculations (including the position of the crossing and the corresponding NAC) can lead to large errors in spectral properties of the molecule. For example, the adiabatic 〈C³Π|DM|X³Σ⁻〉 dipole moment has a steep gradient at around 2 Å which can be expected to be due to the avoiding crossing between C³Π and C′³Π states, therefore the C³Π–X³Σ⁻ electronic band is expected to be sensitive to the quality of its adiabatic description. The diabatic representation can also be sensitive to the quality of the corresponding curves, but to a significantly lesser extend due to their smooth character.

Comparison with the 〈A³Π|μ_x|X³Σ⁻〉, 〈C³Π|μ_x|X³Σ⁻〉, and 〈B³Σ⁻|μ_z|X³Σ⁻〉 transition dipoles provided by Sarka and Nanbu²⁷ shows excellent agreements up to dissociation, with values ({ours, Sarka and Nanbu²⁷}) at the ground state equilibrium geometry R_e(X³Σ⁻) = 1.48 Å of {0.16,0.18} D, {0.333,0.337} D, {1.623,1.633} D, respectively.

3.3 Nuclear motion calculations

Duo⁶³ is a general purpose variational (open access‡) code that solves the rovibronic Schrödinger equation for diatomics while allowing an arbitrary number of couplings between various electronic states including spin-spin, spin–orbit, spin-rotation, and rotational Born–Oppenheimer breakdown curves. It is assumed one has solved the Schrödinger equation for the electronic motion a priori in order to obtain PECs, SOCs, EAMCs, (T)DMCs etc. for the electronic states in question. These curves can be supplied to the program as either a grid of ab initio points, or in an analytical form. After solving the Schrödinger equation for the nuclear motion Duo obtains eigenstates and energies for the good quantum numbers J (total angular momentum), and τ (parity); other quantum numbers are assigned on the basis of the largest coefficient in the basis set. The eigenfunctions are used to compute transition line strengths and Einstein A coefficients in order to obtain a complete spectroscopy for the system in question. A detailed methodology of Duo is given by Duo.

3.4 The ab initio SO spectrum

Using the diabatic spectroscopic model we produce an ab initio rovibronic spectrum of SO for the system involving the lowest 11 singlet and triplet electronic states of SO covering the wavelength range up to 147 nm. The ‘active’ couplings within the spectrosocpic model used for cross-section calculations include 23 SOCs, 23 (transition) DMCs, and 14 EAMCs in-line with the couplings shown in Fig. 2–6. The vibrational sinc-DVR basis set was defined for a grid of 701 internuclear geometries in the range 0.6–6.0 Å. We select 58, 58, 49, 11, 31, 41, 27, 27, 14, 20, and 36 vibrational wavefunctions for the X³Σ⁻, a¹Δ, b¹Σ⁺, A³Π, B³Σ⁻, c¹Σ⁻, A′′³Σ⁺, A′³Δ, C³Π, d¹Π, and e¹Π, respectively, to form the contracted vibronic basis. In total 15 364 624 Einstein A coefficients between 119 600 bound rovibronic states were computed with a maximum total rotational quantum number J_max = 180 and used to simulate rovibronic absorption spectra at a given temperature using the program ExoCross.⁶⁴

Fig. 11 shows an absorption rovibronic spectrum of SO computed at 5000 K with all bands plotted using different colours, both electric dipole allowed and forbidden. We model the spectrum at a high temperature for a visual aid since at this temperature there is a good separation between different electronic bands. The grey shaded region in Fig. 11 marks the total SO bound-bound absorption at 5000 K, which is mostly traced by the strongest bands with the exception for the region between 12 000–17 000 cm⁻¹. Hence, weaker bands, e.g. ones that break dipole selection rules, have negligible contribution to the total SO opacity and will be less important for low resolution studies such as in astrophysical observations.

Fig. 11 Dipole allowed and forbidden components of the ab initio absorption spectrum simulated with the diabatic model at 5000 K connecting the X³Σ⁻, a¹Δ, A′³Δ, and A′′³Σ⁺ states. The C′³Π continuum is also plotted in gold. The absorption lines are modelled using Gaussian profiles with a HWHM of 0.6 cm⁻¹ in the bound cases, and 300 cm⁻¹ for the continuum band.

The non-bound diabatic states such as C′³Π and (3)¹Π are excluded from the bound-bound spectra simulations. Our tests show that the effect of the unbound states on the bound-bound spectra is negligible, and vice versa, the continuum spectra are negligibly affected by the bound electronic states and therefore can be treated separately.

The continuum spectra of the non-bound diabatic C′³Π state is computed using the method described in Pezzella et al.⁶⁵ and is shown in Fig. 11, plotted in gold and overlaying the bound-bound spectrum to demonstrate its contribution to the total SO opacity. For the continuum state a larger basis set of 5000 wavefunctions was used. The structure energetically above (below) the ‘dip’ at 41 200 cm⁻¹ is due to absorption to unbound C′³Π states above (below) the S(¹D) +O(³P) dissociation. The X³Σ⁻ → C′³Π continuum band continues to 100 000 cm⁻¹, peaking at ∼78 000 cm⁻¹ which corresponds to the most vertical transitions from states localised around the minima of X³Σ⁻.

We note that the dipole-forbidden bands in Fig. 11 are not computed using quadrupole or magnetic dipole moments, which have very weak intensities, but rather intensities are ‘stolen’ from other transitions. This intensity stealing propagates through the mixture of electronic wave-functions via couplings such as SOCs and EAMCs. For example, the spin-forbidden c¹Σ⁻–X³Σ⁻ band occurs due to the overlap between the c¹Σ⁻ wavefunction both with e¹Π and d¹Π wavefunctions through the EAM couplings, and then with X³Σ⁻ through a secondary mixing via 〈e¹Π|SO_X|X³Σ⁻〉 and 〈d¹Π|SO_X|X³Σ⁻〉 to produce a direct dipole moment, which dominates over the weaker magnetic and quadrupole moment mechanisms.

We compute absolute intensities for every rovibronic transitions between the lowest 11 diabatic singlet and triplet states of SO covering the entire spectroscopic range up to 147 nm, where NACs are treated. We note that the only other study with similar coverage into the UV on SO is from the theoretical work by Sarka and Nanbu²⁷ who compute cross sections for 190–300 nm. However, our spectroscopic model is both more complete and phase consistent (phases carefully reconstructed, see Section 3.2.2), whereas Sarka and Nanbu²⁷ do not provide any phases.

3.5 Experimental coverage of the ab initio SO spectrum

Currently within the literature a small fraction of the SO spectrum has been measured experimentally covering only the X³Σ⁻, a¹Δ, b¹Σ⁺, A³Π, and B³Σ⁻ states. Fig. 12 reviews the spectroscopic coverage of SO from 24 experimental sources from the literature. Fig. 11 shows our model to supplement the experimental data over the whole spectral range. In particular, we cover the SO spectrum above 40 000 cm⁻¹ and 12 000–16 000 cm⁻¹ where no measurements have been taken for any electronic state. We also plot the available Hitran²⁸ SO line list containing data on the first three electronic states X³Σ⁻, a¹Δ, and b¹Σ⁺. Our ab initio model is able to extend the Hitran coverage up to dissociation at 40 000 cm⁻¹.

Fig. 12 Coverage of experimental measurements for 24 sources illustrated by horizontal bars covering the spectral regions measured, where the named works^6,18,66,67 include some spectral data, mostly with relative intensities: (a) 14 sources^{10,14,15,30,68–78} cover X³Σ⁻ → X³Σ⁻, a¹Δ → a¹Δ, b¹Σ⁺ → b¹Σ⁺ for 0–125 cm⁻¹; (b) 3 experimental sources^17,19,20 cover the A³Π → X³Σ⁻ and B³Σ⁻ → X³Σ⁻ bands for 38 000–39 800 cm⁻¹; (c) 2 experimental sources^21,23 cover the X³Σ⁻ → X³Σ⁻ and a¹Δ → a¹Δ bands for 1040–1125 cm⁻¹.

The aim of a future work is to refine our ab initio SO model to the experimental transition frequencies from these sources and to produce an empirically accurate line list for SO.

4 Effect of diabatisation on the computed spectra

In theory, the adiabatic and diabatic representations should provide identical results, provided that the corresponding NACs and DCs are included. However, due to the computational cost of computing NACs through proper ab initio methods, it is not common practice to include NACs in adiabatic models. Without definition of the NAC, non-Born–Oppenheimer interactions are effectively removed.

In this section we analyse the importance of the non-adiabatic couplings between C³Π and C′³Π as well as between e¹Π and (3)¹Π for computing the (absorption) spectra of SO. To this end we consider our diabatic spectroscopic model to be complete and use it to compute a reference absorption spectrum of SO. This spectrum is then used to compare with an ‘adiabatic’ spectrum of SO computed using the non-adiabatic curves without NACs.

As detailed in Section 3.1, the diabatisation of the potential energy curves induces non-zero off-diagonal elements via corresponding 2 × 2 potential energy matrices in the diabatic representation. These off-diagonal coupling elements will be referred to as diabatic couplings, DCs, which are characterised by cusp shaped curves centered at the avoided crossing. As before, in the diabatic spectrum simulations the unbound states are excluded. Fig. 13 illustrates the vibrational energy levels of the diabatic C³Π, fully bound below its dissociation limit of 50 700 cm⁻¹.

Fig. 13 The computed Duo diabatic vibronic energies of the C³Π superimposed upon the C³Π and C′³Π PECs. The secondary and tertiary ‘bumps’ at R = 2.42, 3.30 Å are due to an avoided crossing and numerical noise, respectively.

Fig. 14 illustrates the importance of the non-adiabatic effects when modeling the spectra around avoiding crossings for the X³Σ⁻ → C³Π and a¹Δ → e′Π bound-bound absorption bands (panel a), and the X³Σ⁻ → C′³Π continuum absorption band (panel b). The adiabatic spectra were computed with the NACs excluded and compared to the diabatic spectra with the non-adiabatic effects fully encountered. Each spectra has been modelled at a temperature of 5000 K – such that hot bands are populated, aiding our comparisons below – with Gaussian profiles of a 0.6 cm⁻¹ half-width-at-half-maximum (HWHM) for the bound-bound spectra and a HWHM of 300 cm⁻¹ for continuum bands.

Fig. 14 A comparison between the X³Σ⁻ → C³Π, a¹Δ → e¹Π, and X³Σ⁻ → C′³Π band spectra computed with an adiabatic model with no defined NAC and a diabatic model. These bands are dipole allowed and are expected to be observable, we see great differences between the spectra at the dissociation, highlighting the importance of diabatisation. Each spectra has been modelled at a temperature of 5000 K with Guassian line profiles of a 0.6 cm⁻¹ HWHM.

Great differences between the bound-bound spectra in panel (a) of Fig. 14 are seen towards both the high and low energy regions. In the high energy region the adiabatic spectral bands terminate abruptly at the avoided crossings whereas the diabatic bands continue to the diabatically correlated dissociation asymptotes S(¹D) + O(³P) & S(¹D) + O(¹D) (see Fig. 1). The diabatisation extends these bands by at least a few thousand wavenumbers because of the availability of higher rovibrational states in the deeper diabatic potential wells. For purely bound-bound calculations, the adiabatically computed bands have lower intensities compared to the diabatic spectrum which can be attributed to the increased repulsive character of the adiabatic PECs on the right hand side of the crossing points present. Due to the tunneling through the potential barriers, the adiabatic wavefunctions ‘leak’ to the continuum region thus resulting in reduction of the intensity of their bound absorption spectra. The most interesting feature from Fig. 14 is the extension of the a¹Δ → e¹Π band beyond the stronger X³Σ⁻ → C³Π band at E/hc > 50 000 cm⁻¹. Although being relatively weak, this band is not covered by stronger bands and therefore may be observable in the ∼0.18–0.2 μm region, a result only predicted when using a full non-adiabatic theoretical treatment.

The low intensity regions are very sensitive to changes in the ab initio model and will be also affected by the changes in the shape of couplings between the adiabatic and diabatic representations. The hump-like structure in the C³Π–X³Σ⁻ band at around 10 000 cm⁻¹ is absent in the adiabatic spectrum because of the unavailability of vibrational states above the avoided crossing. For example, the brightest transitions within this hump for the C³Π–X³Σ⁻ band connect the ν = 13 state which is energetically above the avoided crossing in the adiabatic PEC.

We note that these regions negligibly contribute to the total SO opacity and so are not important for the SO model, but will be important for other systems where non-adiabatic effects occur in the spectroscopically important regions.

Panel (b) of Fig. 14 presents a similar analysis for the continuum X³Σ⁻ → C′³Π band, which would include an additional bound structure towards longer wavelengths if the NAC is not included in the adiabatic model since the C′³Π PEC in this representation is bound. However, the adiabatic X³Σ⁻ → C³Π bound feature is orders of magnitude weaker than the continuum bands presented here and an analysis on the change of character of bound-bound absorption bands with diabatisation is already provided above. The X³Σ⁻ → C′³Π continuum bands for transitions to unbound C′³Π states above the S(¹D) +O(³P) dissociation converge between both representations, since the non-adiabatic effects are far away from the peak at ∼78 000 cm⁻¹ corresponding to vertical transitions from the electronic ground state. However, if the avoided crossing occurred vertically above the X³Σ⁻ minima, we would expect non-adiabtaic effects to have a greater contribution to the continuum cross sections.

From the comparison above, we show that neglecting NACs within an adiabatic model can lead to drastic differences in the physics gleamed from the computed spectra.

5 Conclusions

In this work, we use multireference methods of electronic structure theory combined with a diabatisation procedure to compute a fully diabatic model for the transient diatomic molecule sulphur monoxide. The model includes 23 spin–orbit, 23 (transition) dipole moment, and 14 electronic angular momentum curves for the X³Σ⁻, a¹Δ, b¹Σ⁺, A³Π, B³Σ⁻, c¹Σ⁻, A′′³Σ⁺, A′³Δ, C³Π, d¹Π, and e¹Π electronic states of SO and were produced ab initio via CASSCF and MRCI methods using aug-cc-pV5Z basis sets. These curves were then used to compute the nuclear motion via solving the fully-coupled Schrödinger equation with the Duo program. A further two electronic states (C′³Π and (3)¹Π) were computed along with their couplings, which are essential to forming the diabatic representation. The diabatisation procedure we present is a computationally low cost method to reconstruct the non-adiabatic couplings a priori to nuclear motion calculations. To assess the importance of non-adiabatic effects for the spectroscopy of SO, we compare spectra computed in the diabatic and adiabatic representations without definition of NACs. The most notable difference is the absence of the UV spectrum above ∼50 000 cm⁻¹ because of the illusionary predissociation from the adiabatic PECs. We also saw the adiabatically computed bound absorption bands to have lower intensities than the diabatic counterparts. It is therefore important to treat NACs for systems where these non-adiabatic interactions occur in spectroscopically important regions since they have drastic effects on spectral structure.

All coupling curves of SO are defined with self-consistent relative phases, which is crucial for spectral calculations.³¹ Therefore our spectroscopic model of SO provides a comprehensive and extensive theoretical baseline, which is the first fully reproducible spectroscopic description of SO longer than 147 nm. Since the existing spectroscopic data on SO only covers X³Σ⁻, a¹Δ, b¹Σ⁺, A³Π, and B³Σ⁻, our ab initio model can be used as a benchmark for future rovibronic methods and calculations.

The topic of our next work will be to build a semi-empirical line list as part of the ExoMol project^79,80 for SO through the refinement of our ab initio model to experimental transition data, where we expect to reduce the shift in line positions relative to experiment. The final SO line list will have applications primarily in the atmospheric modelling of exoplanets^81–83 and cool stars. Further applications of this empirical SO line list will be in shock zone modelling,^11,12,84 SO lasing systems,^7,8 and spectroscopy of Venus^85,86 and Io.^87,88

Data availability

The ab initio curves of SO obtained in this study are phase consistent and are provided as part of the supplementary material to this paper along with our spectroscopic model in a form of a Duo input file. The line list computed with Duo can be directly used with the ExoCross program to simulate the spectral properties of SO. We also provide a Julia source code used for the diabatisation procedure is provided as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Appendix

A.1 The Geometric average of DCs

We show here the simple geometric average used to combine the DCs arising from our NACs modelled with a Lorentzian and a Laplacian. We have from eqn (16) of An and Baeck⁴⁸ the following
where in the second line the matrix elements given in eqn (4) are substituted. Comparing the last two lines, one can recover eqn (8), the mixing angle that corresponds to the geometrically averaged DC is then
which corresponds to eqn (15) of An and Baeck.⁴⁸

A.2 Diabatisation procedure

Non-adiabatic interactions cause steep gradients around the region of the avoided crossings in both PECs and coupling curves. To diabatise these couplings one requires the operation of the unitary transformation matrix given in eqn (1). For some coupling, ζ, and considering two adiabatically interacting states |1〉 and |2〉, the exact transformation for non-diagonal couplings 〈1/2|ζ|k〉, where k ≠ 1, 2, from the adiabatic (‘a’) to the diabatic (‘d’) representation can be found as follows
However, diagonal couplings also require the off-diagonal counterparts, and transform similarly to the potential energies as,

The computational procedure for performing the diabatisation is straightforward, and is based on the heuristic that the diabatisation should result in continuous potential energy curves that are twice differentiable due to the properties of the wave function and derivatives. For a function represented by a discrete grid of points V_i at geometries R_i, gradients are calculated via finite differences and the task of satisfying this criterion is approximated by minimizing the sum of second derivatives of the function.

To achieve this we optimize the parameters of an arbitrary NAC function, ϕ(R;{p}). The parameters {p} are typically: a central geometry R_c, and a characteristic width ω (see e.g.eqn (5) and (6)). The NAC function, in turn, parameterises the transformation from the adiabatic potential energy operators V^a_i to the diabatic potential energy operators V^d_i at each geometry i, as described by eqn (4).

The optimization itself is easy to achieve using any commonly available optimization libraries. In the present work we use the Julia programming language and the open-source Optim library's Optim.minimizer function to minimize the following loss function using the Nelder–Mead method:^89,90

where the terms of the summation are simply central finite difference second derivatives of the diabatic potential energy operator.

The integration of the NAC function to obtain the mixing angle β_i at each geometry is obtained by adaptive Gauss-Kronrod quadrature using the QuadGK.jl library.

Empirically we find that an initial guess for the characteristic width ω with the correct order of magnitude suffices for the optimizer to converge. However, the procedure is somewhat more sensitive to the initial guess for the central geometry R_c as a result we provide a convenience function that attempts to detect the central geometry by searching for the largest absolute value of the second order derivatives of the adiabatic operators. The source code for the diabatisation procedure is provided as part of the ESI.†

Acknowledgements

This work was supported by UK STFC under grant ST/R000476/1. This work made use of the STFC DiRAC HPC facility supported by BIS National E-infrastructure capital grant ST/J005673/1 and STFC grants ST/H008586/1 and ST/K00333X/1. We thank the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme through Advance Grant number 883830.

References

1. A. G. Gaydon, Spectroscopy and combustion theory, Chapman and Hall, Ltd., London, 2nd edn, 1948.
2. M. Wang, Encyclopedia of Energy, Elsevier, New York, 2004, pp. 771–789.
3. G. D. Thurston, International Encyclopedia of Public Health, Academic Press, Oxford, 2nd edn, 2017, pp. 367–377.
4. X. Pan, Encyclopedia of Environmental Health, Elsevier, Oxford, 2nd edn, 2019, pp. 823–829.
5. X. Pan, Encyclopedia of Environmental Health, Elsevier, Burlington, 2011, pp. 290–296.
6. J. B. Burkholder, E. R. Lovejoy, P. D. Hammer, C. J. Howard and M. Mizushima, J. Mol. Spectrosc., 1987, 124, 379–392.
7. H. C. Miller, K. Yamasaki, J. E. Smedley and S. R. Leone, Chem. Phys. Lett., 1991, 181, 250–254.
8. B. C. Stuart, S. M. Cameron and H. T. Powell, Chem. Phys. Lett., 1992, 191, 273–278.
9. F. O. Clark and D. R. Johnson, Astrophys. J., 1974, 191, L87–L91.
10. G. Cazzoli, V. Lattanzi, S. Coriani, J. Gauss, C. Codella, A. A. Ramos, J. Cernicharo and C. Puzzarini, Astron. Astrophys., 2017, 605, A20.
11. L. Chernin and C. Masson, Astrophys. J., 1993, 403, L21–L24.
12. M. Y. Amin, M. S. Elnawawy and M. A. Elshalaby, Astrophys. Space Sci., 1991, 185, 277–294.
13. T. Klaus, S. P. Belov, A. H. Saleck, G. Winnewisser and E. Herbst, J. Mol. Spectrosc., 1994, 168, 235–247.
14. S. Yamamoto, Chem. Phys. Lett., 1993, 212, 113–117.
15. M. A. Martin-Drumel, F. Hindle, G. Mouret, A. Cuisset and J. Cernicharo, Astron. J., 2015, 799, 115.
16. R. Colin, Can. J. Phys., 1968, 46, 1539–1546.
17. R. Colin, J. Chem. Soc., Faraday Trans., 1982, 78, 1139–1147.
18. K. D. Setzer, E. H. Fink and D. A. Ramsay, J. Mol. Spectrosc., 1999, 198, 163–174.
19. M. A. A. Clyne and P. H. Tennyson, J. Chem. Soc., Faraday Trans., 1986, 82, 1315–1325.
20. B. C. Stuart, S. M. Cameron and H. T. Powell, J. Phys. Chem., 1994, 98, 11499–11511.
21. H. Kanamori, J. E. Butler, K. Kawaguchi, C. Yamada and E. Hirota, J. Mol. Spectrosc., 1985, 113, 262–268.
22. M. Wong, T. Amano and P. Bernath, J. Chem. Phys., 1982, 77, 2211–2213.
23. H. Kanamori, E. Tiemann and E. Hirota, J. Chem. Phys., 1988, 89, 621–624.
24. P. F. Bernath, R. Johnson and J. Liévin, J. Quant. Spectrosc. Radiat. Transfer, 2021, 272, 107772.
25. P. F. Bernath, R. Johnson and J. Liévin, J. Quant. Spectrosc. Radiat. Transfer, 2022, 108317.
26. L. Yu and W. Bian, J. Comput. Chem., 2011, 32, 1577–1588.
27. K. Sarka and S. Nanbu, J. Phys. Chem. A, 2019, 123, 3697–3702.
28. L. S. Rothman, I. E. Gordon, Y. Babikov, A. Barbe, D. C. Benner, P. F. Bernath, M. Birk, L. Bizzocchi, V. Boudon, L. R. Brown, A. Campargue, K. Chance, E. A. Cohen, L. H. Coudert, V. M. Devi, B. J. Drouin, A. Fayt, J.-M. Flaud, R. R. Gamache, J. J. Harrison, J.-M. Hartmann, C. Hill, J. T. Hodges, D. Jacquemart, A. Jolly, J. Lamouroux, R. J. Le Roy, G. Li, D. A. Long, O. M. Lyulin, C. J. Mackie, S. T. Massie, S. Mikhailenko, H. S. P. Müller, O. V. Naumenko, A. V. Nikitin, J. Orphal, V. Perevalov, A. Perrin, E. R. Polovtseva, C. Richard, M. A. H. Smith, E. Starikova, K. Sung, S. Tashkun, J. Tennyson, G. C. Toon, V. G. Tyuterev and G. Wagner, J. Quant. Spectrosc. Radiat. Transfer, 2013, 130, 4–50.
29. A. Kramida, Y. Ralchenko, J. Reader and NIST ASD Team, Atoms, 2020, 8, 56.
30. C. P. Endres, S. Schlemmer, P. Schilke, J. Stutzki and H. S. P. Müller, J. Mol. Spectrosc., 2016, 327, 95–104.
31. A. T. Patrascu, C. Hill, J. Tennyson and S. N. Yurchenko, J. Chem. Phys., 2014, 141, 144312.
32. M. Born and R. Oppenheimer, Ann. Phys., 1927, 389, 457–484.
33. C. A. Mead and D. G. Truhlar, J. Chem. Phys., 1982, 77, 6090–6098.
34. A. W. Jasper, B. K. Kendrick, C. A. Mead and D. G. Truhlar, Non-Born–Oppenheimer Chemistry: Potential Surfaces, Couplings, and Dynamics, World Scientific, 2004, pp. 329–391.
35. T. Karman, M. Besemer, A. van der Avoird and G. C. Groenenboom, J. Chem. Phys., 2018, 148, 094105.
36. Y. Shu, Z. Varga, A. G. Sampaio de Oliveira-Filho and D. G. Truhlar, J. Chem. Theory Comput., 2021, 17, 1106–1116.
37. J. B. Delos, Rev. Mod. Phys., 1981, 53, 287–357.
38. H. Nakamura and D. G. Truhlar, J. Chem. Phys., 2003, 118, 6816–6829.
39. T. Karman, A. van der Avoird and G. C. Groenenboom, J. Chem. Phys., 2016, 144, 121101.
40. Q. Ma, J. Kłos, M. H. Alexander, A. van der Avoird and P. J. Dagdigian, J. Chem. Phys., 2014, 141, 174309.
41. J. Kłos, Q. Ma, M. H. Alexander and P. J. Dagdigian, J. Chem. Phys., 2017, 146, 114301.
42. T. de Jongh, T. Karman, S. N. Vogels, M. Besemer, J. Onvlee, A. G. Suits, J. O. F. Thompson, G. C. Groenenboom, A. van der Avoird and S. Y. T. van de Meerakker, J. Chem. Phys., 2017, 147, 013918.
43. H. J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, W. Györffy, D. Kats, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, G. Jansen, C. Köppl, Y. Liu, A. W. Lloyd, R. A. Mata, A. J. May, S. J. McNicholas, W. Meyer, M. E. Mura, A. Nicklass, D. P. O’Neill, P. Palmieri, D. Peng, K. Pflüger, R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson and M. Wang, MOLPRO, version 2015.1, a package of ab initio programs, 2015, https://www.molpro.net.
44. T. H. Dunning, J. Chem. Phys., 1989, 90, 1007–1023.
45. D. E. Woon and T. H. Dunning, J. Chem. Phys., 1993, 98, 1358–1371.
46. M. Born and K. Huang, Dynamical Theory of Crystal Lattices, Clarendon Press, 1988.
47. D. Simah, B. Hartke and H.-J. Werner, J. Chem. Phys., 1999, 111, 4523–4534.
48. H. An and K. K. Baeck, J. Chem. Phys., 2015, 143, 194102.
49. K. Baeck and H. An, J. Chem. Phys., 2017, 146, 064107.
50. H. Werner and W. Meyer, J. Chem. Phys., 1981, 74, 5802–5807.
51. A. J. C. Varandas, J. Chem. Phys., 2009, 131, 124128.
52. T. Mondal and A. J. C. Varandas, J. Chem. Phys., 2011, 135, 174304.
53. A. J. C. Varandas, J. Chem. Phys., 2011, 135, 119902.
54. A. Varandas, Chem. Phys. Lett., 2009, 471, 315–321.
55. C.-P. Liu, N. L. Elliott, C. M. Western, Y.-P. Lee and R. Colin, J. Mol. Spectrosc., 2006, 238, 213–223.
56. A. C. Borin and F. R. Ornellas, Chem. Phys. Lett., 2000, 322, 149–156.
57. C. P. Archer, J. M. F. Elks and C. M. Western, J. Chem. Phys., 2000, 112, 6293–6300.
58. A. C. Borin and F. R. Ornellas, Chem. Phys., 1999, 247, 351–364.
59. K. P. Huber and G. Herzberg, Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules, Van Nostrand Reinhold Company, New York, 1979.
60. K. Norwood and C. Ng, Chem. Phys. Lett., 1989, 156, 145–150.
61. H. M. Wang, X. S. Tang, S. K. Zhou, W. J. Zhang and Y. N. Chu, Chem. Phys. Lett., 2005, 407, 78–82.
62. B. Rosen, Spectroscopic Data Relative to Diatomic Molecules, Elsevier Science, 2013.
63. S. N. Yurchenko, L. Lodi, J. Tennyson and A. V. Stolyarov, Comput. Phys. Commun., 2016, 202, 262–275.
64. S. N. Yurchenko, A. F. Al-Refaie and J. Tennyson, Astron. Astrophys., 2018, 614, A131.
65. M. Pezzella, S. N. Yurchenko and J. Tennyson, Phys. Chem. Chem. Phys., 2021, 23, 16390–16400.
66. Y. N. Chu, H. M. Wang, J. Q. Li, P. Cheng and D. Z. Cao, Chem. Phys. Lett., 2002, 366, 147–152.
67. H. M. Wang, X. S. Tang, H. Y. Han, J. Q. Li, S. K. Zhou, W. J. Zhang and Y. A. Chu, Chin. J. Chem. Phys., 2005, 18, 670–674.
68. E. Kim and S. Yamamoto, J. Mol. Spectrosc., 2003, 219, 296–304.
69. F. X. Powell and D. R. Lide, J. Chem. Phys., 1964, 41, 1413–1419.
70. M. Winnewisser, W. Gordy, K. Sastry and R. L. Cook, J. Chem. Phys., 1964, 41, 1687–1691.
71. E. Tiemann, J. Mol. Spectrosc., 1974, 51, 316–320.
72. W. W. Clark and F. C. Delucia, J. Mol. Spectrosc., 1976, 60, 332–342.
73. E. Tiemann, J. Mol. Spectrosc., 1982, 91, 60–71.
74. Y. Endo, H. Kanamori and E. Hirota, Chem. Phys. Lett., 1987, 141, 129–132.
75. G. Cazzoli, L. Cludi, G. Cotti, C. D. Esposti and L. Dore, J. Mol. Spectrosc., 1994, 167, 468–471.
76. T. Klaus, A. H. Saleck, S. P. Belov, G. Winnewisser, Y. Hirahara, M. Hayashi, E. Kagi and K. Kawaguchi, J. Mol. Spectrosc., 1996, 180, 197–206.
77. M. Bogey, S. Civis, B. Delcroix, C. Demuynck, A. F. Krupnov, J. Quiguer, M. Y. Tretyakov and A. Walters, J. Mol. Spectrosc., 1997, 182, 85–97.
78. T. Klaus, S. P. Belov and G. Winnewisser, J. Mol. Spectrosc., 1997, 186, 416–421.
79. J. Tennyson and S. N. Yurchenko, Mon. Not. R. Astron. Soc., 2012, 425, 21–33.
80. J. Tennyson, S. N. Yurchenko, A. F. Al-Refaie, V. H. Clark, K. L. Chubb, E. K. Conway, A. Dewan, M. N. Gorman, C. Hill, A. Lynas-Gray, T. Mellor, L. K. McKemmish, A. Owens, O. L. Polyansky, M. Semenov, W. Somogyi, G. Tinetti, A. Upadhyay, I. Waldmann, Y. Wang, S. Wright and O. P. Yurchenko, J. Quant. Spectrosc. Radiat. Transfer, 2020, 255, 107228.
81. M. Y. Zolotov and B. J. Fegley, Icarus, 1998, 132, 431–434.
82. R. Hobbs, P. B. Rimmer, O. Shorttle and N. Madhusudhan, Mon. Not. R. Astron. Soc., 2021, 506, 3186–3204.
83. V. A. Krasnopolsky, Icarus, 2012, 218, 230–246.
84. R. Bachiller, Annu. Rev. Astron. Astrophys., 1996, 34, 111–154.
85. C. Y. Na, L. W. Esposito and T. E. Skinner, J. Geophys. Res.: Atmos., 1990, 95, 7485–7491.
86. D. A. Belyaev, F. Montmessin, J.-L. Bertaux, A. Mahieux, A. A. Fedorova, O. I. Korablev, E. Marcq, Y. L. Yung and X. Zhang, Icarus, 2012, 217, 740–751.
87. E. Lellouch, Icarus, 1996, 124, 1–21.
88. I. de Pater, H. Roe, J. R. Graham, D. F. Strobel and P. Bernath, Icarus, 2002, 156, 296–301.
89. J. Bezanson, A. Edelman, S. Karpinski and V. B. Shah, SIAM Rev., 2017, 59, 65–98.
90. P. K. Mogensen, J. M. White, A. N. Riseth, T. Holy, M. Lubin, C. Stocker, A. Noack, A. Levitt, C. Ortner, B. Johnson, D. Lin, K. Carlsson, Y. Yu, C. Rackauckas, J. Grawitter, A. Williams, B. Kuhn, B. Legat, J. Regier Cossio, R. Rock, T. R. Covert, B. Pasquier, T. Arakaki, A. Stukalov, A. Clausen, A. Strouwen and B. Deonovic, JuliaNLSolvers/Optim.jl: v1.7.2, 2022, https://zenodo.org/record/7019119.

---

Footnotes

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d2cp03051a

‡ https://github.com/Exomol/Duo)github.com/Exomol.

---