ExoMol line lists -- LVI: The SO line list, MARVEL analysis of experimental transition data and refinement of the spectroscopic model

A semi-empirical IR/Vis line list, SOLIS, for the sulphur monoxide molecule $^{32}$S$^{16}$O is presented. SOLIS includes accurate empirical rovibrational energy levels, uncertainties, lifetimes, quantum number assignments, and transition probabilities in the form of Einstein $A$ coefficients covering the $X\,{}^{3}\Sigma^{-}$, $a\,{}^{1}\Delta^{ }$, $b\,{}^{1}\Sigma^{+}$, $A\,{}^{3}\Pi$, $B\,{}^{3}\Sigma^{-}$, $X\,{}^{\prime\prime3}\Sigma^{+}$, $A\,{}^{\prime 3}\Delta$ and $e\,{}^{1}\Pi$ systems and wavenumber range up to 43303.5 cm$^{-1}$ ($\geq 230.93$ nm) with $J\le 69$. SOLIS has been computed by solving the rovibronic Schr\"{o}dinger equation for diatomics using the general purpose variational code Duo and starting from a published ab initio spectroscopic model of SO (including potential energy curves, coupling curves, (transition) dipole moment curves) which is refined to experimental data. To this end, a database of 50106 experimental transitions, 48972 being non-redundant, has been compiled through the analysis of 29 experimental sources, and a self-consistent network of 8558 rovibronic energy levels for the $X$, $a$, $b$, $A$, $B$, and $C$ electronic states has been generated with the MARVEL algorithm covering rotational and vibrational quantum numbers $J \leq 69$ and $v \leq 30$ and energies up to 52350.40 cm$^{-1}$. No observed transitions connect to the $ B\,{}^{3}\Sigma^{-} (v = 0)$ state which is required to model perturbations correctly, so we leave fitting the $B\,{}^3\Sigma^-$ and $C\,{}^3\Pi$ state UV model to a future project. The SO line list is available at ExoMol from www.exomol.com.


INTRODUCTION
Sulphur monoxide ( 32 S 16 O) is an often transient diatomic molecule whose spectral properties are important for a wide range of environments.SO was first detected by radio astronomy in the interstellar medium by means of rotational spectroscopy and was the first 3 Σ ground state molecule to be detected in outer space (Gottlieb & Ball 1973).Additionally, SO was the first molecule observed using pure rotational transitions using microwave spectroscopy within its excited electronic states (Saito 1970).Since its first detection, SO has been observed in many astronomical environments, including the interstellar medium (Gottlieb et al. 1978), molecular clouds (Codella & Muders 1997;Blake et al. 1987), and planetary and lunar atmospheres (Lellouch 1996;de Pater et al. 2002;Na et al. 1990;Belyaev et al. 2012).Numerous studies propose sulphur-bearing molecules, including SO, as constituents of volcanic planetary atmospheres (Zolotov & Fegley 1998;Hobbs et al. 2021;Krasnopolsky 2012).It also plays a role in many solar-system atmospheres, including that of Jupiter's moon Io (Lellouch 1996;de Pater et al. 2002) and of Venus (Na et al. 1990;Belyaev et al. 2012), as well as Earth's own atmosphere (Burkholder et al. 1987).The recent detection of SO 2 in at atmosphere of WASP-39b (Rustamkulov et al. 2023) suggested SO as part of the photochemical production paths of SO 2 (Tsai et al. 2023).SO also has importance for: (1) experimental applications in spectroscopy, such as UV lasing (Miller et al. 1991;Stuart et al. 1992); (2) ★ The corresponding author: j.tennyson@ucl.ac.uk astrophysical applications, such as shock modelling (Bachiller 1996;Chernin & Masson 1993;Amin et al. 1991), magnetic field measurements in molecular clouds through observation of SO Zeemen splitting (Clark & Johnson 1974;Cazzoli et al. 2017), and studies into planetary formation mechanisms after its first detection within a protoplanetary disk (Pacheco-Vazquez et al. 2015), which supports sulphur-carrier based dust-formation channels (Laas & Caselli 2019;Vidal et al. 2017); (3) telluric applications, such as environmental studies of acid rain, since it is an intermediate in combustion reactions and has great chemical involvement with with N 2 and O 2 (Gaydon 1948;Burkholder et al. 1987).
The rovibrational structure of SO's spectrum for the  3 Σ −  1 Δ,  1 Σ + ,  3 Π, and  3 Σ − electronic states at low vibrational excitation's has been studied by numerous works; we provide a full analysis of the experimental coverage on SO below.SO's electronic transitions were first reported by Martin (1932), and has since been subject to pure rotational (Klaus et al. 1994;Yamamoto 1993), electronic (Colin 1968(Colin , 1982;;Setzer et al. 1999), and ro-vibrational (Kanamori et al. 1985;Wong et al. 1982;Kanamori et al. 1988) spectroscopic studies.The pure rotational transitions within several vibrational states in the ground  3 Σ − electronic state have been measured in the terahertz (Martin-Drumel et al. 2015;Klaus et al. 1997), far infrared (Clark & Delucia 1976;Cazzoli et al. 1994), and microwave regions (Yamamoto 1993;Lovas et al. 1992;Bogey et al. 1997).For most of these spectra, SO was studied in non-local thermodynamic equilibrium (non-LTE) conditions, and so only relative intensities are available at best.Currently, there are no published absolute inten-sity measurements for SO.However, measurements of state lifetimes provide information on Einstein A coefficients and hence transition dipole moments (Tennyson et al. 2016).Experimental lifetimes for the  1 Σ + ,  3 Π and  3 Σ − states have been measured (Elks & Western 1999;Smith 1969;Yamasaki et al. 2003), and provide a valuable benchmark for our intensity calculations (see Section 5.1).
In this study we present the largest compilation of experimental transition data and derived self-consistent empirical rovibrational energy levels for 32 S 16 O to date.The derived energy levels where obtained through use of the MARVEL (Measured-Active-Rotational-Vibrational-Energy-Levels) spectroscopic network algorithm, to which we format the 32 S 16 O data based on the MARVEL format (Furtenbacher et al. 2007).We then refine our ab initio spectroscopic model (Brady et al. 2022) to our determined empirical energies to produce a hot semi-empirical line list SOLIS for 32 S 16 O as part of the ExoMol project (Tennyson & Yurchenko 2012;Tennyson et al. 2020).The SOLIS line list supplements existing spectroscopic line list data for SO which are limited in coverage.For example, spectroscopic databases CDMS (Endres et al. 2016) and NIST (Kramida et al. 2020) databases contain data covering the microwave region only.HITRAN (Gordon et al. 2022) considers relatively low vibrational excitations for transitions between electronic states  3 Σ − ,  1 Δ and  1 Σ + only.We compare the SOLIS line list to the existing spectra data in Section 5.5.

The MARVEL Procedure
The critical evaluation of experimental transition data and formation of a self-consistent set of rovibronic energy levels is done through the MARVEL procedure (Furtenbacher et al. 2007;Furtenbacher & Császár 2012a) which is built on the concept of spectroscopic networks (SN) (Furtenbacher & Császár 2012b;Császár & Furtenbacher 2011).Through a weighted linear least squares protocol, MARVEL inverts the information contained within transition data to form a set of associated energy levels and uncertainties.Self consistency within the energy levels is then achieved through an iterative re-weighting algorithm which adjusts (increases) the uncertainties in the line positions to an optimised uncertainty  opt until they agree with the rest of the network.
The inverted MARVEL energy levels form nodes of a SN, which are linked by transitions, to which the validation of experimental information can be done on all data simultaneously using elements of network theory.The final energy level uncertainties in the SN are obtained through combining the optimal MARVEL uncertainties of all transitions connecting a given energy level.This study used a new implementation, MARVEL 4, which uses a bootstrap method to determine uncertainties in the empirical energy levels it determines (Tennyson et al. 2023).We used 100 iterations with the bootstrap method to determine the uncertainties in our empirical MARVEL energy levels.

Quantum Numbers
We assign the rovibronic energy levels of 32 S 16 O using the vibrational and rotational quantum numbers  and , respectively, rotationless parity  (/  ) and, in line with Hund's case-(a) coupling scheme and sublevels denoted by the fine structure  2+1 .For  ≥ Λ + , the fine structure,  2+1 , is defined for triplet electronic states via where singlet states have no spin projection (i.e. 2 ) and the total angular momentum excluding electronic and nuclear spin is labelled .For linear molecules such as diatomics, levels with parity +(−1)  and −(−1)  are labelled  and  levels, respectively, and their relation to the ± parities are given in Table 1.We thus assign every experimental rovibrational transition using the , , /  ,  quantum numbers and standard spectroscopic notation for electronic states.
For nuclear motion calculations we use the quantum numbers Λ, Σ, and Ω to assign electronic states, which are the projection of orbital, spin, and total angular momentum on the bond axis, respectively, in additional to the state labels  3 Σ − ,  3 Π,  1 Σ + etc.For the  3 Σ − and  3 Σ − states the spin-parity sub-levels in increasing energy order are whereas for the regular  3 Π state there is lambda-doubling and we have and for the inverted  3 Π state the sublevels increase in energy with decreasing Ω Rigorous electric-dipole selection rules hold here, and can be summarised as

Outline
Table 2 summarises the experimental transition data included within our MARVEL analysis where each study is conveniently labelled with a tag including the first two digits of the year of publication and letters of the names of the first three authors in the form 'YYAaBbCc'. Table 2 includes the spectral coverage of each study, the associated quantum number coverage of their assignments, and the mean uncertainty of their results.We compiled a total of 50 106 transitions, of which 49 613 are non-redundant, from 29 experimental sources covering the  3 Σ − ,  1 Δ,  1 Σ + ,  3 Π,  3 Σ − and  3 Π electronic states of SO for rovibrational excitation's  ≤ 69,  ≤ 30.

General Comments
A crucial limitation of the experimental data set for SO (Table 2) is in the vibrational state coverage of the lower electronic states.
Transitions to/within states beyond the third vibrational excitation for  3 Σ − ,  1 Δ,  1 Σ + , and  3 Π are severely lacking.Some vibrational transition data involving states beyond  > 3 are available, but inclusion of these within our MARVEL analysis often led to fragmented SNs.Our literature review found that no transitions have been measured associated with the vibrational ground state of the  3 Σ − state, which makes it difficult to constrain its PEC minima during refinement of the spectroscopic model.The impact of this missing data is then amplified since the  3 Σ − and  3 Π energies exhibit many mutual perturbations because of their overlapping potentials and strong coupling.To correctly model the perturbations one requires accurate positioning of the potentials corresponding to the resonating states relative to each other, which is made difficult because of the absence of data connecting to the  3 Σ − ( = 0) state.

Source Specific Comments
(a) A significant problem we faced during data analysis is that several literature sources did not provide obvious uncertainties on their line measurements, which is important for their validation within the MARVEL protocol (Section 2.1).We thus had to estimate their uncertainties through combination difference (CD) tests to other data in our database with known uncertainties where possible.These sources include Bouchoux et al. (1971a); Bouchoux & Marchand (1972);Colin (1969Colin ( , 1982) ) and Stuart et al. (1994).Initial uncertainties for these sources were assumed to be 0.05 cm −1 and manually increased with successive MARVEL runs until the data gave satisfactory CD relations with other sources.In the case of blended lines, their uncertainties were increased by a factor of √ 2 relative to the non-blended data which often resulted in their validation.As a result of this the source uncertainties were estimated to be 0.02 cm −1 (Stuart et al. 1994), 0.4 cm −1 (Bouchoux et al. 1971b(Bouchoux et al. , 1972)), 0.05 cm −1 (Colin 1969) and 0.2 cm −1 (Colin 1982).
(b) Another issue with the experimental transition data is the significant proportion of blended lines that are reported, such in the measurements by Colin (1969Colin ( , 1982)); Clyne & Tennyson (1986); Burkholder et al. (1987); Kanamori et al. (1988) and Bogey et al. (1997).To account for potential inaccuracy in their assignments the blended lines were given a lower weight in our SN model (see comment (a) above).
(d) 15MaHiMo (Martin-Drumel et al. 2015) provide much data on SO isotopologues with determined isotopically invariant parameters as well as other various constants for the lowest 7 vibrational states.
(f) 88KaTiHi (Kanamori et al. 1988) provide 60 SO ( 1 Δ,  = 3− 4, 4−5) transitions.If high vibrational data for SO becomes available these would be a prime source for inter-vibrational transition data to supplement our MARVEL dataset.
(g) A discrepancy in the CDMS data for the  3 Σ − ,  = 1 state of SO was found.The lowest  = 1 state energy had to be shifted by 26.4559 cm −1 .Furthermore, a shift of 6.478 cm −1 was found in the  1 Δ,  = 0 CDMS data of SO, where the source of error may come from use of a high uncertainty  = 0 − 0 band centre due to Bielefeld et al. (1984).The relative energies between CDMS levels within the same vibrational levels are unaffected, hence the transition wavenumbers are correct, but were corrected before being used in analyses involving MARVEL.
(h) Heays et al. (2022) perform high-resolution FUV Fouriertransform photoabsorption spectroscopy and provide the only published UV transition data covering the  3 Π← 3 Σ − and  3 Σ − ( = 4 . . .30) ←  3 Σ − bands.Because of the large overlap and spin-orbit coupling (SOC) between the  3 Σ − and  3 Π states, many perturbations are present within the experimental data which appear to be assigned accurately.They also provide transition data for the isotopologues 33 S 16 O and 36 S 16 O.
(i) Setzer et al. (1999) provide 74  1 Σ + −  3 Σ − magnetic dipole transitions, which have the same selection rules as for electric dipole transitions except from the parity changing rule.We do not include these in our MARVEL network.
(j) We chose to omit the 540  3 Π ( ′ = 2)− 3 Σ − ( ′′ = 0) transitions measured by Colin (1969) for two reasons: (1) they produced many conflicts with the more comprehensive and more accurate data by Heays et al. (2022); (2) use of MARVELised energies generated including these data for refining our spectroscopic model proved very difficult; we prescribed abnormalities in the energy structure to be due to the poor data which did not occur when using the equivalent MARVELised data from Heays et al. (2022).

MARVELisation of the Experimental Transition Data
The 32 S 16 O spectroscopic network was built through input of 50 106 rovibronic transitions into MARVEL from the 29 data sources outlined in Table 2. 546 transitions were invalidated since their optimised uncertainties did not satisfy the validation condition  opt −  exp < 0.05 cm −1 .Invalidation of transition data can be due to multiple reasons, errors in their quantum number assignment, in their measurement, in the digitization of their scanned data Tables (especially in old papers), and simply because they are not self-consistent with the rest of the network.The latter is the most common cause of inval-Table 2. The experimental data sources included in the final MARVEL analysis and their spectroscopic coverage.TAG denotes the identifier used to label the data sources throughout this paper, V/T describes the number of validated or included (V) data using the MARVEL procedure described in Section 2.1 relative to the total number of provided transitions (T), and the final column cross-references source specific comments (Comm) in Section 3.3.idation but since they usually connect few levels they are invalidated if the aforementioned reasons are not the cause.The invalidated transitions are removed from the MARVEL network but are kept in the MARVEL input file with a negative wavenumber transition frequency.We note the 590 transitions were excluded not because of invalidation through the MARVEL procedure and are detailed in comments (i) and (j) in Section 3.3.The majority of transitions that were invalidated are for lines connecting  > 3 (69%) because of the lack of inter-vibrational data energetically above  > 3 which, if included, resulted in the fragmentation of our central SN and the invalidation of otherwise seemingly reliable data sources.There are a lack of measurements of rotational transitions within these higher vibrational states.It was found that keeping data for  ≤ 3 produced the largest set of self consistent energy levels and hence SN.For  ≤ 3 the experimental source that provided the most invalidated transitions is by Colin (1969) (44%) who measured the only  3 Π − →  3 Σ −  = 2 − → 0 band transitions, which is important for the refinement of the  3 Π potential energy curve (PEC).Colin (1969) measured 32 S 16 O in emission by means of flash photolysis of sulphur bearing gases using a medium resolution quartz spectrograph.They provide no direct uncertainty on their line positions, but provide an uncertainty for their  3 Σ − − →  3 Σ − bandheads of ±1 cm −1 obtained in their absorption study.If one uses this value as a metric for their line position uncertainties, then it is to be expected that data coming from Colin (1969) should be treated with caution, consequently leading to much of their data being invalidated.The majority of the remaining invalidated  ≤ 3 data comes from Colin (1982) and Stuart et al. (1994) (21% and 22%, respectively) As a result of the critical evaluation of the experimental transition data, we invert and provide optimised uncertainties for 8558 rovibrational energy levels for 32 S 16 O which forms a fully self-consistent SN. Figure 1 plots the MARVEL energies versus the rotational quantum number , where a large gap in the ∼15 000-37 000 cm −1 region exists, corresponding to missing highly excited vibrational data in the  3 Σ − ,  1 Δ, and  1 Σ + states and any experimental coverage of the intermediate electronic states  1 Σ − ,  ′ 3 Δ, and  ′′ 3 Σ + .For the higher vibrational energy levels of each electronic state there are also gaps in the rotational structure.The experimental transition frequencies collected as part of this work are provided in the Supporting Information to this paper in the MARVEL format alongside the resulting MARVEL empirical energy levels.

THE SPECTROSCOPIC MODEL
To produce the final semi-empirical line list for 32 S 16 O we use the ab initio spectroscopic model presented in our recent work (Brady et al. 2022) as a theoretical baseline for refinement to our MARVELised energy levels.Section 4.1 overviews the details of the ab initio model and Section 4.2 details the method used to refine it.

The Ab initio Spectrosopic Model
The ab initio curves that make up our SO model (Brady et al. 2022) were computed using internally contracted multireference configuration interaction (ic-MRCI) level of theory and aug-cc-pV5Z basis sets.The active space and state averaging were chosen to have 12 active electrons with occupied (8,3,3,0) and closed (4,1,1,0) orbitals  Figure 1.Our generated MARVEL energies plotted against the angular momentum quantum number  for the  3 Σ − , 1 Δ,  1 Σ + ,  3 Π,  3 Σ − , and  3 Π states.The vertical structure within each electronic state corresponds to the different vibrational levels.The size of the plot markers are directly proportional to the number of combination differences to that level in the SN.
under C 2 symmetry.The model covers 13 electronic states: and (3) 1 Π, which range up to 66800 cm −1 , vastly beyond the scope of experimental coverage of the molecule.The ab initio model includes 13 potential energy curves, 29 dipole and transition dipole moment curves, 25 spin-orbit curves, and 18 electronic angular momentum curves on a grid of points over bond lengths 1-3 Å, where a diabatic representation was built by removing the avoided crossings between the spatially degenerate pairs  3 Π- ′ 3 Π and  1 Π-(3) 1 Π through a property-based diabatisation method.

Refinement
We refined our model to the IR/Vis region by fitting to the  3 Σ − ,  1 Δ,  1 Σ + , and  3 Π energies only.We include the minimum number of states and couplings required such that our computed energies for these states agree with the MARVEL energies.Our refined model includes PECs and couplings connecting the 2, where potentials for states labelled with a ' †' superscript are not refined; these are included solely for their couplings within our model, but their dipoles are kept ab initio.These couplings, despite not being the dominant contributions to the energies of the  3 Σ − ,  1 Δ,  1 Σ + and  3 Π states, will redistribute intensities according to the intensity stealing mechanism (see Sec.5 and the Appendix).Firstly, we found that including the  1 Σ − and  1 Π states within the model had negligible effect on the energies of the states of interest, and so they were omitted.Secondly, the inclusion of the  3 Π state proved problematic and only worsened the fit.Because the  3 Σ − and  3 Π PECs have large overlap and strong couplings between each other, their MARVEL energies include many perturbations because of resonances with each other, and also due to multiple avoided crossings with upper electronic states.Furthermore, due to the lack of important experimental data and consequent lack of proper constraint on the  3 Σ − PEC minimum, efforts to include the UV region within our fit proved too difficult to do satisfactorily.We therefore removed the  3 Π state entirely from our model, but found including the  3 Σ − , now without resonances with  3 Π, improved the fit of the  3 Σ − ,  1 Δ,  1 Σ + , and  3 Π energies without being a major contribution to the accuracy of our model.Our initial efforts to fit the UV model constrained the  3 Σ − PEC enough such that expectation values of its couplings to other states in our model were sensible.Furthermore, the current need for the IR/Vis SO line list means we leave work on the UV model to a future study when the appropriate data becomes available.
Refinement of the ab initio model is facilitated through Duo, a general purpose variational (open access 1 ) code that solves the rovibronic Schrödinger equation for diatomic molecules.A description of the methodologies used in Duo is given by Yurchenko et al. (2016).
The refinement process goes as follows: (1) represent PECs, SOCs, EAMCs, DMCs, and other empirically fitted couplings such as rotational Born-Oppenheimer breakdown, spin-spin, and spin-rotational curves with analytical forms; (2) compute energy levels using Duo through solving the rovibronic Schrödinger equation with curves defined; (3) fit parameters of the analytical functions such that the computed energy levels agree with the MARVEL energies.
Before refinement, the MARVEL energy level quantum number assignments need to be converted to the Duo quantum numbers , , , Λ, Σ, and Ω (see Section 2.2).Next, the -th MARVEL energy was given a weight equal to where is its optimised uncertainty and  CD  is the number of combination differences/frequency of occurrence within the transition database.This weight is used in the Duo fit meaning energy levels with large uncertainty and a low number of combination differences will have less effect on the fit.Next, a running number must be defined to enumerate the global order of the energy levels calculated by Duo.These should agree with MARVEL's energy ordering for lower  and  numbers, but for higher energy states where Duo's calculated energies deviate significantly from the MARVEL data, the ordering of states between the two can differ.The running number we employ increases by 1 per vibrational excitation and by 100 per each electronic state, for example the running numbers for  3 Σ − ( = 1) and  1 Δ( = 2) are 1 and 102, respectively.This produced a sensible enumeration that effectively separated the energy levels and provided correct assignments of the calculated energy levels.

Potential energy, spin-orbit, electronic angular momentum curves
We represent all PECs using the Extended Morse Oscillator (EMO) function (Lee et al. 1999) which has the form where  e =  e −  e is the dissociation energy,  e is the potential minimum and  e is the asymptote,   are the expansion coefficients,  e is the equilibrium bond length, and   () is the so called Šurkus variable (Šurkus et al. 1984) given by with  as an integer parameter to allow for a better convergence at large bond lengths values.The and  ′′ 3 Σ + states dissociate to the same asymptote S( 3 P) + O( 3 P), which we initially set at 5.429 eV as reported by Huber & Herzberg (1979) and then floated during our fits, which converged to a nearby value of 5.42895 eV.The  3 Σ − state adiabatically correlates to S( 1 D) + O( 3 P) which we set to a value of 6.5731 eV (53015.86cm −1 ) as determined through atomic energies from the NIST atomic database.We shifted the PECs to the  3 Σ − minimum such that  e ( 3 Σ − ) = 0 cm −1 .Since only the  3 Σ − ,  1 Δ,  1 Σ + ,  3 Π,  3 Σ − and  3 Π have experimental transition data to refine their PECs to, the ab initio PECs for the  ′′ 3 Σ + , and  ′ 3 Δ were fitted to the EMO function given in Eq. ( 6) using 10 expansion parameters which ensured accurate representation of their shape as given by ab initio calculations.Once we fitted the EMO functions, we could then tune their dissociation assymptotes to the 5.429 eV limit (Huber & Herzberg 1979).The  e value for the  ′′ 3 Σ + state was fixed to 30 692 cm −1 as provided by Norwood & Ng (1989).We chose not to tune the  1 Π ab initio PEC because it has a strong influence on the computed  3 Σ − ,  1 Δ,  1 Σ + , and  3 Π energies and negatively effects our refinement when altering its potential.This negative effect is due to: (1) The shape of the PEC would be destroyed in tuning the  e and  e values; (2) the  1 Π state has been diabatised (Brady et al. 2022), so tuning its PEC would change the avoided crossing morphology and hence a new diabatisation of the spectroscopic model would be required; without experimental data covering  1 Π we chose to keep the diabatised ab initio potential values.
In our refinement the ab initio  3 Σ − PEC by Sarka & Nanbu (2019) was used instead of the PEC from our recent ab initio work (Brady et al. 2022).The latter did not employ sulphur specific diffuse functions, which led to the underestimation of the S( 3 P) + O( 3 P), S( 1 D) + O( 3 P) and S( 1 D) + O( 1 D) dissociation asymptotes and effected the  3 Σ − PEC the worst.The adiabatic character of the  3 Σ − PEC was modelled by diagonalising a 2 × 2 matrix of diabatic potentials and the corresponding diabatic coupling which gives the adiabatic potential through solution of the quadratic equation for this system (Brady et al. 2022).
The ab initio SOCs and EAMCs were morphed from the grid representation to a Šurkus-like expansion (Prajapat et al. 2017;Yurchenko et al. 2018a) given by where   are the expansion coefficients,  ∞ is usually taken as zero in-order to allow the expansion to not diverge towards  − → ∞, and  is a damped displacement coordinate given by where  2 and  4 are damping constants.

Empirical Rotational Born-Oppenheimer Breakdown, Spin-Spin, Spin-Rotational Curves
We fitted the phenomenological spin-spin (SS) couplings and the empirical spin-rotation (SR) couplings of the triplet  3 Σ − ,  3 Π, and  3 Σ − states to account for additional Ω-splitting and to allow for additional variation of , respectively, not described by the ab initio model (Yurchenko et al. 2016;Kato 1993).We also fited rotational Born-Oppenheimer breakdown (BOB) curves for all but the  3 Σ − state to correct for the electron un-coupling to the nuclear motion producing an additional  2 dependence in the residuals of the rovibronic energies and which can be thought of as a correction to the position-dependent rotational mass.Some SS, SR, and BOB couplings are modelled using Eq. ( 9) and some using a Šurkus polynomial expansion given by This greatly enhanced the accuracy of the finalised spectroscopic model.

Accuracy of the Refined Model
Figure 3 illustrates the Observed minus Calculated (Obs.-Calc.)energy residuals as a function of the rotational quantum number , and provides a metric on the accuracy of our model to reproduce our MARVELised energies.Most of the highly scattered energy levels have no combination differences with other sources, and so are effectively removed from the fit by setting their weight to 10 −6 .We fit 100% of 512  3 Σ − ( ≤ 69) energy levels with a total root-meansquare (rms) error of 3.13 × 10 −3 cm −1 , 99% of 244  1 Δ ( ≤ 52) energy levels with a rms error of 1.08 × 10 −3 cm −1 , 95% of 206  1 Σ + ( ≤ 64) energy levels with a total rms error of 0.27 cm −1 , and 78% of 1262  3 Π ( ≤ 34) energy levels with a rms error of 0.24 cm −1 .The above rms errors were calculated after filtering outliers from our dataset which heavily influenced the rms, such as the scattered data of  1 Σ + ( = 3, 4), the  3 Π( = 2,  ≥ 35) states, and a single data point of  1 Δ.Elks & Western (1999) provide rotational constants for the  3 Π ( = 4 − 13) states fitted to their 1+1 resonance enhanced multiphoton ionisation spectra which we used to compute low  ( ≤ 5) energies via Pgopher Western (2017).These energies to our knowledge are the only ones covering the highly excited vibronic states of  3 Π and so we use them to constrain the  3 Π potential up to its dissociation.The black points in figure 3 for  ≤ 5 (see label (a)) indicates the vibrational dependence in the residuals of these Pgopher levels, which we manage to fit all within ∼4.5 cm −1 .Doing this allowed for a more physical description of the effective position dependent correction to the rotational mass for the  3 Π by constraining the potential gradient, and ultimately led to higher accuracy in the associated computed lifetimes (see Section 5.1).Generally, the  1 Σ + data are of high quality and are reproduced by our model to within ∼ 10 −4 − 10 −3 cm −1 with the exception of high scatter within the MARVELised  1 Σ + ( = 3, 4) energies.This scatter can be seen in Fig. 3 around the label (b).One striking feature in Fig. 3 is of the  3 Π ( = 2) Obs.-Calc.residuals, which, despite many attempts to model correctly, are poorly recovered in our model.We see a smooth, but rapid increase in the MARVEL energies with  of this band to ∼35 cm −1 at  = 51.We postulate this is because a dark state pushes these energies upwards, some candidates being the  1 Σ − ,  ′′ 3 Σ + , and  ′ 3 Δ which cross through the  3 Π potential, however all attempts to correctly position their PECs relative to the  3 Π failed to reproduce this behaviour.No published data on the crossing states exists, but it is entirely possible that a correct description of these electronic potentials could resolve this issue, and some empirical data could be used to constrain their curves better than blindly varying their positions.
Some residual -dependence can be seen for the  3 Σ − and  1 Σ + states, where the former is due to -dependent Ω-splitting and parity splitting within the |Ω| = 1 levels, and the latter is due to vibrational dependence in the effective rotational centre not being fully accounted for.However, the residuals to the MARVEL energies are all ≤10 −2 cm −1 which means the model should extrapolate well to higher .
One major problem faced during the refinement was with the spinorbit splitting of the  3 Π energies, where experiment (Colin 1969(Colin , 1982) ) predicts regular Ω energy ordering, whereas multiple ab initio calculations reveal the  3 Π SO z A 3 Π SOC to have a negative phase, suggesting irregular Ω energy ordering.Analysis by Colin (1969) shows the Λ-doubling (LD) to be ∼1.2 cm −1 in the lowest energy Ω state with slight dependence on , small doubling in the middle coupling which varies with  ( + 1) and zero splitting for the highest component state.To choose whether we adopt the irregular Ω energy ordering suggested by ab initio calculations or the experimental assignment with a change in phase of the  3 Π SO z A 3 Π SOC we studied the /  parity splitting of the  3 Π energies since this would confirm what Λ-doubling matrix elements to adopt.We saw that the slightly -dependant ∼1.2 cm −1 parity splitting in the lowest energy state could only be resolved via the following Λ-doubling with the Λ-doubling constant being consistent with the Brown & Merer (1979) convention  LD   =  LD +  LD +  LD for a regular Ω energy ordering, i.e. the splitting is between states of ΔΣ = ±2 which is only possible for a triplet state if Ω = 0.If one adopted an irregular assignment scheme with Ω = 2 for the lowest energy state, then one could not correctly model the -dependence of the splitting since the element in Eq. ( 11) would be zero.We thus adopted the experimental assignment and changed the sign of our  3 Π SO z A 3 Π SOC which should not break the phase consistency of the model since it is a diagonal coupling.

Dipole Moment Curves
We use the accurate ab initio ground state dipole moment function from Bernath et al. (2022) who calculate with an ic-MRCI+Q level of theory including the Davidson corrections, scalar relativistic contributions using the exact 2-component (X2C) relativistic Hamiltonian, and aug-cc-pCV6Z-X2C basis sets.All other DMCs are computed at a level of theory described in Section 4.1.Within nuclear motion and intensity calculations, these dipoles are originally represented as a grid of ab initio points on the Duo defined grid, however one sees a flattening of both the IR  3 Σ − - 3 Σ − band spectrum and its variation of TDMC with vibrational excitation.The source of this nonphysical flattening has been discussed by Medvedev et al. (2015) and Medvedev et al. (2016); Medvedev & Ushakov (2022) who identify numerical noise as the culprit.This noise comes from Duo interpolating the given Molpro dipole grid points onto the Duo defined grid.One can try to increase the precision of their transition moments from double to quadruple precision, but this seldom fixes the problem with any appreciable magnitude.The most effective  2022), computed at an ic-MRCI+Q level of theory with full relativistic corrections using aug-cc-pCV6Z-X2C basis sets, is shown (black crosses) superimposed with our fitted analytical form using Eq. 12 (red line).The residuals to the ab initio DMC of our fit are shown in the top panel (green line).method found is to represent the input dipole moments analytically.We chose to represent our  3 Σ − DMC using the 'irregular DMC' proposed in Medvedev & Ushakov (2022) Where   are Chebyshev polynomials of the first kind,   are summation coefficients to be fit, () is a reduced variable in bond length similar to the damped polynomial coordinate in Eq. ( 9) and is given by, which maps the  ∈ [0, ∞] interval to the  ∈ [−1, 1] reduced interval, and finally (;  2 , ...,  6 ) is an -dependent term parametrically dependent on 5   parameters to be fitted and is given by Our fitted  3 Σ − DMC is illustrated with its residual to the ab initio DMC in Fig. 4.
The irregular DMC has the desirable properties of quickly converging to the correct long-range limit, having enough parameters (13) to ensure an accurate description of the full range in bond length with minimal local oscillations, and provide a straight Normal Intensity Distribution Law (NIDL) (Medvedev & Ushakov 2022;Medvedev 2012;Medvedev et al. 2015).
The PEC of the  3 Π state implies that predissociative and continuum states should exist for the region above dissociation.To this end, these states have been removed from the line list through checking the character of the wavefunctions at the 'right' simulation border  max in our Duo model where unbound states tend to oscillate at  → ∞ with a non-zero density around  max (Yurchenko et al. 2022).The line list we present therefore only contains bound to bound transitions only.
The calculated energies in the .statesfile are 'MARVELised' which involves replacing them with the MARVEL ones.For levels that are not covered by the MARVEL SN, the predicted shift method of Bowesman et al. (2021) was used to MARVELise them.Predicted shifts work by fitting the Obs.-Calc.trends as functions of  for each 'state',  and Ω energy band to then interpolate gaps within the MARVEL network or extrapolating to higher .
The SOLIS line list is available in the ExoMol database (https://www.exomol.com/) in the form of a States (.states) and Transition (.trans) files, with extracts shown in Tables 3 and 4 respectively.Uncertainties for the energy levels where either taken directly as the MARVEL ones where available, or otherwise computed using the following empirical formulae where  is the energy uncertainty for a given state and Δ, Δ, Δ are state dependent parameters given in Table 5. Δ were found by taking twice the standard deviation of the total Obs.-Calc. of each electronic state (see Fig. 3) after outliers where removed by selecting states outside of this two standard-deviation threshold, where the standard deviation was computed again.

Intensity Scaling: Dipoles and Lifetimes
There are only a few recorded experimental values for electric dipole moments, lifetimes, and no direct intensity measurements which can be used to constrain our ab initio dipoles.Lifetimes are useful to constrain dipole moments via the relation where l is the dominant ro-vibronic state contributing to the lifetime of the level |⟩,  = 0, ±1 denotes a tensorial dipole component.So a scaling in lifetime τ =  would correspond to an approximate scaling in dipole moment to the dominant lower state of 1/ √ .Previous Stark measurements (Powell & Lide 1964;Lovas et al. 1992) have determined the ground state dipole to be  0  = 1.55(2)D (Powell & Lide 1964) and 1.52(2) D (Lovas et al. 1992), slightly smaller than our computed vibrational transition moment of |  3 Σ − ,  = 0  0  3 Σ − ,  = 0 | = 1.588D. We scale our  3 Σ − dipole to the value of 1.535 D averaged from the two Stark measurements, which we note is the dipole moment adopted by CDMS  (Endres et al. 2016) where the same averaging was done.Wildt et al. (1983) measure the radiative lifetime of the  1 Σ + ( = 0) state through time-resolved measurements of the  1 Σ + − →  3 Σ − emission band and provide a lifetime of  = 6.8 ± 0.4 ms.To achieve this lifetime, we scale our  1 Σ +  0  1 Σ + dipole by a factor of 0.7401.Saito (1970) determines the  1 Δ dipole moment to be 1.336 ± 0.045 D through Stark measurements, larger than our computed transition moment |  1 Δ,  = 0  0  1 Δ,  = 0 | = 1.184D. We scale our  1 Δ dipole by a factor of 1.1282 to reproduce the measured transition moment.
Radiative lifetimes of the  3 Π state for  ′ = 0−13 were measured by Elks & Western (1999) by laser induced fluorescence and Clyne & Liddy (1982) for  ′ = 0 − 6. Figure 7 shows the experimentally determined lifetimes as a function of  ′ with the theoretically predicted values by Borin & Ornellas (2000) and Fulscher et al. (1992) and those computed by our model superimposed in red.Since Elks & Western (1999) quote their lifetimes to a lower uncertainty and for a large vibrational coverage we chose to model their lifetimes.Modelling these lifetimes proved to be difficult, the characteristic sharp drop in lifetime from  ′ = 0 and  ′ = 1 was very sensitive to multiple factors: (1) the position of the  3 Π  ±1  3 Σ − dipole relative to the respective PECs; (2) the crossing point of the dipole with zero; (3) the local gradient of the dipole around the zero crossing point.Initial attempts to reproduce the experimental lifetimes were made using our ab initio dipole, various ab initio dipoles from the literature (Fulscher et al. 1992;Sarka & Nanbu 2019;Feng & Zhu 2019), and the empirical dipole from Elks & Western (1999) in our model all failed to produce lifetimes that agreed with experiment.
Firstly, we were able to reproduce the lifetimes of Elks & Western (1999) using a linear dipole function.Albeit being nonphysical, it provided important constraints on the short range position of and its gradient around the equilibrium geometry.We then combined it with the MRCI-F12+Q/aug-cc-pV(5+d)Z dipole computed by Sarka & Nanbu (2019) at larger values of  into a single smooth curve.Despite a slight change in the shape of the ab initio dipole, Fig. 7 shows that our semi-empirically fitted  3 Π  ±1  3 Σ − dipole is much closer to the ab initio curve than the empirically fitted dipole by Elks & Western (1999).

Partition Function
We compute the molecular partition function (PF) for SO from our semi-empirical line list using where  2 is the second radiation constant, Ẽ is the rovibronic energy term value in wavenumbers,  tot  =  ns (2  + 1) is the total state degeneracy which includes the nuclear weight spin-statistic  ns ( ns =  Sauval & Tatum (1984), Barklem & Collet (2016), CDMS (Endres et al. 2016), and Hitran (Gamache et al. 2017) who compute their PF from the line lists produced by Bernath et al. (2021Bernath et al. ( , 2022)).As the nuclear spin degeneracy is one, no PFs need to be scaled to the physics convention of nuclear statistical weights, which ExoMol uses. Figure 8 shows that all PFs agree for 500 ≲  ≲ 2000 K; for all temperatures our computed PF continues to agree with that of Barklem & Collet (2016), where our computed PF is generally lower than theirs up to 0.1% at 5000 K; the CDMS PF agrees to within 1% of all PFs up to it cutoff at 300 K; Sauval & Tatum (1984) is larger than the rest of the PFs at low temperatures up to 500 K and at 5000 K their PF is 3.4% lower than our computed PF; the Hitran PF begins to deviate from the other PFs from about 2000 K, where at 5000 K it is lower than ours by 17%; this behaviour of Hitran PFs has been noted previously (Wang et al. 2023).

SO opacities
We follow the ExoMolOP procedure of Chubb et al. (2020) and generate molecular opacitites for SO using the SOLIS line list for four exoplanetary atmosphere retrieval codes ARCiS (Min et al. 2020), TauREx (Al-Refaie et al. 2021), NEMESIS (Irwin et al. 2008) and petitRADTRANS (Molliére et al. 2019) on an extensive grid of temperatures and pressures.The opacities are provided as part of the SO ExoMol data set at www.exomol.com.

Simulated spectra
Program ExoCross (Yurchenko et al. 2018b) was used to simulate rovibronic absorption spectra as a function of temperature using SO-LIS.Figure 5 illustrates the dipole allowed and forbidden electronic bands connecting the  3 Σ − state to  3 Σ − ,  1 Δ,  1 Σ + , and  3 Π which are shown as different colours and the total computed SO opacity is shown in grey.Here we simulate lines with a Gaussian line profile of HWHM 0.6 cm −1 .The forbidden band intensities are stolen through mixing of the electronic wavefunctions through couplings such as SOCs, DMCs, and EAMCs resulting in non-zero dipole matrix elements, which we note provides a stronger mechanism here than their corresponding magnetic dipole or electric quadrupole couplings.Figure 6 shows the temperature variation of the simulated total We see the intensity deviation is greatest in region  around 18 000-35 000 cm −1 where the  3 Σ − ← 3 Σ − band begins to dominate opacity.SO opacity which has a strong effect on the UV/Vis cross-sections.
The greatest temperature variation can be seen in the 18 000-35 000 cm −1 region (V) where the  3 Σ − − →  3 Σ − band begins to dominate opacity.Here we simulate lines with a Gaussian line profile of HWHM 0.6 cm −1 .It is clear the IR/NIR spectrum is largely unaffected by the increase of temperature except from the expected rotational broadening.Below we comment on the spectral regions marked by I-VI illustrated in Fig. 5.
(III) The ∼7800-13000 cm −1 NIR region is dominated by strong  1 Σ + ←  3 Σ − band absorption for all temperatures, and is almost as strong as the dipole allowed  3 Σ − ←  3 Σ − band spectrum because of large intensity stealing mechanism facilitated through the strong  1 Σ +    3 Σ − SOC, see the Appendix.
(IV) The Vis ∼13000-18200 cm −1 region shows a flat feature due to  3 Π←  3 Σ − band absorption which becomes prominent for temperatures above 3000 K.However, since we omit the  3 Π state in our spectroscopic model we do not compute  3 Π←  ′ 3 Δ and  3 Π←  ′′ 3 Σ + band intensities which we previously predicted to be strong in this region (Brady et al. 2022).
(V) The Vis/UV ∼18200-37700 cm −1 region is largely uncovered at high accuracy by our spectroscopic model since it is dominated by the  3 Σ − ←  3 Σ − and lesser  3 Π←  3 Σ − electronic bands which become major sources of SO opacity for temperatures above 1000 K.We are currently working on the UV SO line list for a future study which will accurately cover this region.However, for lower temperatures the  3 Π←  3 Σ − and  1 Σ + ←  3 Σ − bands become more important which we recover accurately in this study.
(VI) The UV ∼37700-43500 cm −1 region has a strong  3 Π←  3 Σ − band feature which should be observable at all temperatures.rovibrational band for 0-7000 cm −1 .We simulate the spectra using a temperature of 296 K and scale the intensities by the fractional isotopologue abundance of 0.9479 (Gordon et al. 2022).

Comparisons to experimental spectra
There are few recorded experimental spectra of 32 S 16 O with large coverage and almost none with absolute intensity measurements, relative intensities are usually provided (Chu et al. 2002;Burkholder et al. 1987;Setzer et al. 1999;Wang et al. 2005).The only study to our knowledge that provides measured absolute intensities is the recent study by Heays et al. (2022) on the  3 Π− 3 Σ − band, which we compare to (also  3 Σ − − 3 Σ − and  3 Π− 3 Σ − bands which we do not compare to).
The forbidden band intensities we compute here are through the intensity stealing mechanism which works through mixing of electronic state wavefunctions through couplings such as SOCs.We do not compute magnetic dipole intensities, which are much weaker for the bands of interest than the redistributed intensities which we compute.For example, the diagonal  1 Σ +    1 Σ + and  3 Σ −    3 Σ − dipoles produce  1 Σ + ← 3 Σ − band intensities a factor of ∼ 10 − 1 000 times stronger at the band peak than the corresponding magnetic dipole intensities.Therefore, we omit magnetic dipole transitions from our line list.An example of the intensity stealing mechanism is given in Appendix.

Hitran
HITRAN produces empirical SO line lists which have been produced by fitting spectroscopic models to experimentally derived spectroscopic constants, lifetimes, and rotational branching ratios (ratios in perpendicular and parallel transition moments, see below dis- cm −1 .We simulate the spectra using a temperature of 296 K and scale the intensities by the fractional isotopologue abundance of 0.9479 (Gordon et al. 2022).
cussion), hence we will compare to this data since it is the closest comparison of rovibronic intensities to experimental data for the  3 Σ − − 3 Σ − ,  1 Σ + − 3 Σ − , and  1 Δ− 3 Σ − electronic bands.Hitran (Gordon et al. 2022) provides empirical line list data on the first three electronic states of SO  3 Σ − ,  1 Δ, and  1 Σ + to which we compare our theoretical spectra to.The HITRAN intensities for the  3 Σ − ← 3 Σ − band were originally presented by Bernath et al. (2022) and the forbidden  1 Σ + ←  3 Σ − and  1 Δ ←  3 Σ − bands are from Bernath et al. (2021).In both studies, fitted spectroscopic constants from the literature were used to predict line positions, transition moments were obtained using LeRoy's LEVEL program (Le Roy 2017) which assumes the single state approximation, and their line lists were computed using PGOPHER (Western 2017).Bernath et al. (2021) used Ω-representation to allow for the single state approximation in line with LEVEL such that the forbidden band intensity are computed from effective dipoles between single Ω-states, the so-called parallel and perpendicular transition moments, as opposed to the non-approximate intensity stealing mechanism via mixing of electronic wavefunctions through, e.g.SOCs, as we do.Perpendicular and parallel electronic transition moments between the spin-orbit states 0 + − 0 + , 0 + − 1, and 2 − 1 (see appendix) were computed by Bernath et al. (2021) at an iC-MRCI/aug-cc-pCVQZ-DK level of theory and were scaled to the experimentally determined values by Setzer et al. (1999).The HITRAN  3 Σ − ← 3 Σ − intensities were computed using the ab initio ground state expectation dipole moment computed by Bernath et al. (2022) at a ic-MRCI+Q/ACV6Z-X2C/ED+Q level of theory.
In all comparisons below, we scale our computed intensities with the 32 S 16 O isotopologue abundance 0.947926 given by Hitran (Gordon et al. 2022).Figure 9 presents a comparison between our semiempirical  3 Σ − ←  3 Σ − rovibronic spectrum, where we compute stick spectrum using a temperature of 296 K, and the empirical Hitran line list (Bernath et al. 2022).We see good agreement in both the line positions and band structure where band intensities agree up to the fifth hot band at ∼ 4500 cm −1 where we compute higher intensities relative to the Hitran data.The 0 ← 0 band agrees extremely well which can be seen in the sub plot of Fig. 9.The agreement in intensities confirms our methodology since the Hitran dipole was also scaled to the same experimental values discussed in Section 5.1.The discrepancy in intensities towards hotter bands can be attributed to the difference in the DMCs as well as the wavefunctions used to calculate the transition probabilities.
Figure 10 compares our computed  1 Σ + ←  3 Σ − spectrum simulated at a temperature of 296 K to the empirical Hitran (Gordon et al. 2022)  1 Σ + ←  3 Σ − line list (Bernath et al. 2021).For this comparison, we filtered out the magnetic dipole transitions present in the Hitran line list, since we only calculate electric dipole transitions.The selection rules for magnetic dipole branches are the same except they follow the non-parity changing rule.We see that our model supplements the Hitran line list at both the higher and lower wavenumber regions (< 4000 cm −1 and > 12500 cm −1 ) where line positions, band structure and intensities generally show good agreement which is to be expected since the HITRAN  1 Σ + ←  3 Σ − dipoles were scaled to the same experimental values for the  1 Σ + lifetime as we do.For vibrational bands energetically below ∼ 7500 cm −1 the two spectra begin to deviate from each other, where SO-LIS tends to be slightly lower than the Hitran intensities.However, this is below the standard Hitran intensity threshold of 10 −30 cm 2 /molecule, and so these bands are typically of less importance.For intensities above the threshold, we see good agreement in line with the methodologies used.
We also compare to the forbidden  1 Δ ←  3 Σ − band in Fig. 10 to the theoretical Hitran line list at 296 K.The transitions here are all electric dipole in nature.We see the electronic band structure agrees well between the peak band up to the penultimate hot band before the Hitran terminus at ∼11000 cm −1 with the largest intensity deviation being between vibronic bands outside the 5000-8500 cm −1 spectral range.The low energy bands <4500 cm −1 all sit below the 10 −30 cm 2 /molecule intensity threshold, which are less important spectroscopically.Differences in the band intensity are difficult to disentangle since  1 Δ ←  3 Σ − is a dipole forbidden band where intensities are accumulated through 'intensity stealing' mechanism via multiple coupling channels in our model.
The general agreement with the empirical Hitran line list confirms we are not missing any fundamental physics since we both produce similar spectra using a different methodology.Furthermore, Bernath et al. (2021) scale their  1 Σ + ←  3 Σ − effective dipoles using the same lifetime provided by Setzer et al. (1999) as we did, whereas the  1 Δ ←  3 Σ − band have no reliable experimental dipoles or lifetimes to scale the ab initio dipoles to except the indirect  1 Δ- 1 Δ dipole moment presented by (Saito 1970) which only we scale to.Agreement for the latter band then confirm the good quality of the PECs, SOCs, and (T)DMCs used.

A-X Bands
The recent study by Heays et al. (2022) measured, in high resolution, the FUV  3 Π− 3 Σ − band via Fourier-transform spectroscopy up to the ( ′ , ") = (3, 0) band for  ≤ 51; the  3 Σ − − 3 Σ − and  3 Π− 3 Σ − bands were also measured but we do not compare to these here.Heays et al. present an empirical line list where effective Hamiltonian spectroscopic constants were fitted band-by-band to their measured spectrum, providing quantum number assignments and oscillator strengths for each assigned transition.Their coupledband models reproduce the experimentally measured line positions, intensities, and widths to within 5% uncertainty.With this, we converted their line list to the ExoCross format to compute corresponding cross-sections, which we compare to.For all spectra simulations a temperature of T = 360 K and Gaussian line broadening of HWHM 0.3 cm −1 was used.Figure 11 shows the comparison between our computed  3 Π( ′ = 1, 2, 3)− 3 Σ − ( ′′ = 0) band intensities in blue (bottom panel) and the simulated band intensities of Heays et al. (2022) in red (top panel).We see excellent agreement in line positions, intensities, and band structure, where a mirror plot was chosen since overlaying the spectra made it hard to distinguish between them since they agree so well.

CONCLUSION
We present the semi-empirical SOLIS line list for 32 S 16 O constructed starting from the refinement of the ab initio spectroscopic model presented by Brady et al. (2022) to empirically derived energy levels, or MARVEL energy levels.As part of the line list creation, a MAR-VEL analysis of 29 experimental transition sources resulted in a self-consistent set of 8558 rotation-vibration energy levels ( ≤ 69 and  ≤ 3) for the  3 Σ − ,  1 Δ,  1 Σ + ,  3 Π,  3 Σ − ,  3 Π electronic states, where 48 972/50 106 experimental transitions were validated.The SOLIS SO line list supplements existing ExoMol line lists for SO 2 (Underwood et al. 2016a) and SO 3 (Underwood et al. 2016b).
The  3 Σ − state expectation value of the dipole moment operator was fitted to an analytical form and shown to improve the nonphysical flattening of the vibrational transition moment and NIDL compared to using the grid interpolated form of the DMC.This resulted the physical exponential decay of the  3 Σ − − →  3 Σ − IR spectral band.
Comparison of the simulated rovibronic spectrum of SO to experiment/Hitran show good agreements in both positions and intensities.However, inspection of the  = 0 → 0 of the forbidden  1 Σ + − →  3 Σ − band revealed disagreeing P-and R-branch ratios to experiment, where tuning of the spectroscopic model showed no proper inversion of the branch intensities.Analysis of the electronic wavefunctions revealed that the band intensities are dominated by competing diagonal  and  DMCs, which contribute to the parallel transition moment.The weaker perpendicular transition moment was shown to produce intensities of the desired P-and R-branch ratio, but was much weaker than the parallel transition component to the intensities.Analysis on the basis set revealed the branch ratio to be more sensitive when including a larger vibrational basis, but still does not produce the desired branch ratio.
The future work includes extension to the UV region with the  3 Σ − ←  3 Σ − and  3 Π←  3 Σ − electronic bands, and production of photodissociation cross sections and rates.where  is the rotational quantum number,  is the parity, and  represents the full set of quantum numbers |⟩ = |State, , Ω, Λ, , Σ, ⟩.
We see that because of large SOC between  3 Σ − and  1 Σ + , they share sizeable contributions in their final mixed wavefunctions of their corresponding unmixed basis states.The amount of intensity stealing will then distribute itself through subsequent coupling of the dipole operators in the new mixed state basis.To this end, we consider the parallel and perpendicular transition dipole moments (TDMs) which couple Ω = 0 + − 0 and Ω = 0 + − 1 states, respectively.In spherical tensor form they read, considering the  3 Σ − and  1 Σ + states, The experimental measurement and analysis by Setzer et al. (1999) of the  1 Σ + − →  3 Σ − emission band and work by (Bernath et al. 2021) shows the  0 and  ±1 TDMs to be of the same order of magnitude.However, comparing the  1 Σ + − →  3 Σ −  = 0 − 0 band measured by Setzer et al. (1999) to our semi-empirical line list and the Hitran empirical line list show disagreements in the P-and R-branch ratios where experiment predicts the P-branch have the most intensity within the band.Considering the intensity stealing mechanism, cross-examining Eqs.(18,19) with Table 6 reveals that intensity stealing is different for the perpendicular and parallel transitions and therefore the branch ratios for these transitions.This could explain the experimental observation in Setzer et al. (1999) that the P-and R-branch ratios are different.In an attempt to understand this discrepancy, our analysis on the mixed  3 Σ − and  1 Σ + state wavefunctions reveals that the competition between the  1 Σ +    1 Σ + and  3 Σ −    3 Σ − DMCs provides the dominant contributions to  0 , as facilitated through the large  1 Σ + SO z X 3 Σ − spin-orbit coupling, where  1 Σ +    1 Σ + subtracts from  3 Σ −    3 Σ − and can be seen in Figure 12 as the reduction of intensity between the red and blue spectra. 1 , however, has leading contributions from the  3 Π    3 Σ − dipole as facilitated through the  1 Σ + SO x A 3 Π and  3 Σ − SO x A 3 Π couplings.We find that the perpendicular TDM  1 is responsible for the P-R-branch asymmetry ( > ) as seen in the Fourier transform spectroscopy by Setzer et al. (1999) whereas the parallel TDM  0 produces the opposite branch ratios ( < ).However, when considering only a vibrationless expansion for our wavefunctions, that is  = 0, the  1 TDM is very weak producing much lower intensities than the parallel transition moment and as a result the  1 Σ + − →  3 Σ − band has a  <  branch ratio.When considering a larger vibrational basis ( = 20) the P-R-branch ratio becomes more sensitive to changes in DMCs and SOCs.From this analysis we see that the vibrational TDMs have an effect on the rovibrational band intensity distribution in favour of the P-R-branch ratio as seen in Setzer et al. (1999).

Figure 2 .
Figure 2. Potential energy curves of states included within our spectroscopic model of SO.The solid lines correspond to the potentials refined to MARVEL data, dashed lines correspond to states included within our model but have not been refined, where couplings to these states are essential for the accuracy of the  3 Σ − ,  1 Δ,  1 Σ + , and  3 Π energies and band intensities.

Figure 4 .
Figure 4.The ab initio  3 Σ − DMC provided by Bernath et al. (2022), computed at an ic-MRCI+Q level of theory with full relativistic corrections using aug-cc-pCV6Z-X2C basis sets, is shown (black crosses) superimposed with our fitted analytical form using Eq. 12 (red line).The residuals to the ab initio DMC of our fit are shown in the top panel (green line).

Figure 5 .
Figure 5. Dipole allowed and forbidden components of the absorption spectrum simulated with our semi-empirical model at 1000 K connecting  3 Σ − with  3 Σ − ,  1 Δ,  1 Σ + , and  3 Π.Regions of spectral importance are marked with roman numerals and are detailed in the text.The grey shaded region marks the total SO opacity computed with our model at 1000 K.

Figure 6 .
Figure6.The total absorption spectrum of SO simulated with our semiempirical model for different temperatures ranging from 298 K to 5000 K.We see the intensity deviation is greatest in region  around 18 000-35 000 cm −1 where the  3 Σ − ← 3 Σ − band begins to dominate opacity.
Σ − and  1 Σ + wavefunctions in the Ω representation.The column headers with a subscript 'Duo' are the computed Duo-states which have components in the basis states given as rows due to spin-orbit coupling. 1 Σ + − →  3 Σ − BAND,AN EXAMPLE OF INTENSITY CONTRIBUTIONS FROM A FORBIDDEN BANDDipole forbidden transitions can arise through multiple mechanisms, such as through the magnetic dipole moment (MDM), quadrupole moment, and from intensity stealing.Intensity stealing propagates through the mixture of electronic state wave-functions via couplings such as SOCs and EAMCs, where contributions to the forbidden intensities are derived through taking dipole matrix elements in the eigenstates of the diagonalised Hamiltonian constructed from the coupled Λ − Σ basis.To understand these intensity contributions to the  3 Σ − − →  1 Σ + band we performed an analysis of the Duo computed electronic state wavefunctions corresponding to the eigensolutions of the diagonalised Hamiltonian which included SOCs, EAMCs, and DMCs for the full 11 state system described inBrady et al. (2022) plus additional  1 Σ + SO x A 3 Π ,  1 Δ SO x A 3 Π and  1 Σ + SO x C 3 Π couplings.The contributions to the  3 Σ − and  1 Σ + computed wavefunctions in the Ω representation are shown in Table6which gives the expansion coefficients   of the wavefunctions in the eigenbasis of the diagonalised Hamiltonian Ψ

Figure 12 .
Figure 12.Visualisation of the different contributions to the forbidden  3 Σ − − →  1 Σ + band intensities.'Q' refers to the magnetic dipole moment, where the green spectra are due to magnetic dipole transitions, which is orders of magnitude weaker than the intensity stealing mechanism.

Table 1 .
Relation between the /  and ± parities for linear molecules.

Table 4 .
Extract from the transitions file of the line list for SO.   : Einstein- coefficient in s −1 ; ν   : transition wavenumber in cm −1 .

Table 5 .
State dependent parameters (in cm −1 ) of Eq. (14) used to estimated uncertainties for the calculated states of 32 S 16 O where MARVEL uncertainties were not available.
16O) where we use a 1 K temperature step.Figure8compares our computed PF to the PFs of

Table 6 .
The largest expansion coefficients for the  3