Spectroscopy of YO from first principles

We report an ab initio study on the spectroscopy of the open-shell diatomic molecule yttrium oxide, YO. The study considers the six lowest doublet states, $X\,{}^{2}\Sigma^{+}$, $A\,{}^{2}\Pi$, $C\,{}^{2}\Pi$, $A'\,{}^{2}\Delta$, $B\,{}^{2}\Sigma^{+}$, $D\,{}^{2}\Sigma^{+}$ and a few higher-lying quartet states using high levels of electronic structure theory and accurate nuclear motion calculations. The coupled cluster singles, doubles, and perturbative triples, CCSD(T), and multireference configuration interaction (MRCI) methods are employed in conjunction with a relativistic pseudopotential on the yttrium atom and a series of correlation-consistent basis sets ranging in size from triple-$\zeta$ to quintuple-$\zeta$ quality. Core-valence correlation effects are taken into account and complete basis set limit extrapolation is performed for CCSD(T). Spin-orbit coupling is included through the use of both MRCI state-interaction with spin-orbit (SI-SO) approach and four-component relativistic equation-of-motion CCSD calculations. Using the ab initio data for bond lengths ranging from 1.0 to 2.5 A, we compute 6 potential energy, 12 spin-orbit, 8 electronic angular momentum, 6 electric dipole moment and 12 transition dipole moment (4 parallel and 8 perpendicular) curves which provide a complete description of the spectroscopy of the system of six lowest doublet states. The Duo nuclear motion program is used to solve the coupled nuclear motion Schr\"{o}dinger equation for these six electronic states. The spectra of $^{89}$Y$^{16}$O simulated for different temperatures are compared with several available high resolution experimental studies; good agreement is found once minor adjustments are made to the electronic excitation energies.


Introduction
Oxides of transition metals and lanthanides have rich and complex spectra due to the presence of many low-lying excited electronic states. This complexity poses particular challenges for experimental 1 and theoretical 2 studies. The yttrium oxide, YO, is an example of a rare-earth oxide whose electronic structure is very difficult to explore. Yttrium is a relatively abundant rare-earth element both on Earth (the 28th most abundant element 3 ) and in space (the second most abundant rare-earth metal 4 ). As a result, the spectrum of YO has been the subject of many astrophysical observations. In particular, YO has been observed in a variety of spectra of cool stars including R-Cygni 5 , Pi-Gruis 6 , V838 Mon 7,8 , and V4332 Sgr 7 . The spectrum of YO has been extensively used as a probe to study high temperature materials at the focus of a solar furnace [9][10][11] . The A 2 Π 1/2 electronic state YO has a relatively short life time of 33 ns 12 with large diagonal Franck-Condon factors, 13 which makes this molecule well suited for cooling experiments with the potential in quantum information applications. 14 Yttrium oxide is one of the very few molecules that have been laser 0 [15][16][17][18] .
YO is a strongly bound system. The compilation by Gaydon 51 reports its dissociation energy to be 7.0±2 eV, while Ackermann and Rauh 52 recommended a D 0 value of 7.290(87) eV based on mass spectrometric determinations.
A few theoretical investigations of YO are available in the literature. The most comprehensive one was carried out by Langhoff and Bauschlicher 53 who reported the spectroscopic constants for the lowest five doublet, X 2 Σ + , A 2 ∆, A 2 Π, B 2 Σ + , C 2 Π, and fourteen quartet electronic states of YO. The doublets were studied at the multireference single and double excitation configuration interaction (MRCI) level of theory and, in the case of the X 2 Σ + , A 2 ∆, and A 2 Π states, also using the modified coupled-pair functional (MCPF) method. All the quartet states were considered at the CASSCF level, and that with the lowest energy, reportedly 4 Φ, at the MCPF level as well. Zhang et al. 37 have recently reported the CASPT2 spectroscopic constants and excitation energies for a set of lowest doublet states of YO including the D 2 Σ + state in addition to the doublets studied previously by Langhoff and Bauschlicher 53 . In all of the previous theoretical studies, only modest double-ζ 53 or triple-ζ 37 basis sets were employed. RKR curves and some Franck-Condon factors of YO were computed by Sriramachandran and Shanmugavel 54 .
The main objective of the present study is to characterize both the electronic ground state and the plethora of low-lying excited states of YO with high-level ab initio methods, and to accurately describe from first principles the spectroscopy of YO via producing the potential energy curves (PECs) and other data needed to calculate the rovibronic energies and transition probabilities comprising a so-called line list for this molecule. The generation of such line lists is a major object of the ExoMol project. 55 Thus far, ExoMol studies of open-shell transition metal (TM) diatomics have struggled due to difficulties in providing reliable ab initio starting points. 2,[56][57][58] The intrinsic challenge to theory posed by open-shell systems is associated with several types of problems including spin contamination, symmetry breaking in the reference function, strong nondynamical electron correlation effects, avoided crossings between adiabatic potential energy surfaces, etc. (for the discussion, see, e.g., Refs. 59,60). In the openshell TM-containing species, these problems are exacerbated by stronger relativistic effects than those in the molecules made up of relatively light main group elements, and greater number of electronic excited states governing the spectroscopic behaviour of a molecule and hence deserving to be taken into account in a study aimed at accurate description of its spectroscopy. Moreover, the low-lying electronic states of TM species are commonly degenerate or near-degenerate, which complicates their theoretical treatment even more. Multireference methods of quantum chemistry best suited for describing closely spaced electronic states might seem to be the natural choice for studying these systems. However, most routine multireference methods, such as MRCI, are incapable of properly handling dynamical electron correlation and therefore do not provide high accuracy description of TM-containing species commonly featuring strong dynamical correlation effects. Such effects are best treated with single reference coupled cluster (CC) theory known for its capability to predict highly accurate properties even for molecules with mild to moderate MR character. Unfortunately, the higher likelihood of severe multireference character in the ground and/or low-lying electronic excited states of open-shell TM-containing species makes their treatment by single reference methods very problematic, if possible at all. Particularly this is true for the studies aimed at a description of the molecular potential energy surfaces over a wide range of geometries. It is therefore not surprising that the high-level coupled cluster studies on the open-shell TM-containing species, where a few excited states are treated on an equal footing with the ground state, are very uncommon and only deal with near-equilibrium regions of these states (see, e.g., Refs. [61][62][63][64]. Such a study on a manifold of electronic excited states of a TM-containing diatomic molecule over a wider bond length range has not been reported so far.
It is thus clear that none of routine methods of modern quantum chemistry are entirely satisfactory in all respects for accurately describing from first principles the spectroscopy of open-shell TMcontaining species. Nevertheless, one can try to solve this challenging task via the so-called composite approach by which the desired set of molecular properties is obtained using multiple methods of different nature and sophistication rather than a single method.
In this paper, we have examined efficiency of such an approach taking the example of YO. The PECs for the six lowest doublet electronic states of this molecule, X 2 Σ + , A 2 ∆, A 2 Π, B 2 Σ + , C 2 Π, D 2 Σ + , were obtained from the extensive high-level coupled cluster calculations addressing core-valence correlation and basis set convergence issues, whereas the spin-orbit curves (SOCs), electronic angular momentum curves (EAMCs), electric dipole moment curves (DMCs), and transition dipole moment curves (TDMCs) were obtained at the MRCI level of theory. These curves, with some simple adjustment of the minimum energies of the PECs, are used to solve the coupled nuclear-motion Schrödinger equation with the program DUO. 65 The spectroscopic model and ab initio curves are provided as part of the supplementary material. Our open source code DUO can be accessed via http://exomol.com/software/.

Ab initio calculations
Multireference single and double excitation configuration interaction, MRCI, 66-68 and underlying complete active space self-consistent field, CASSCF, 69,70 calculations were carried out using a relativistic energy-consistent 28-electron core pseudopotential (PP) accompanied with the aug-cc-pVTZ-PP 71 basis set on the Y atom, and the aug-cc-pVTZ 72 all-electron basis set on the O atom (this combination of sets is hereafter referred to as aVTZ). To obtain a consistent MRCI data set in the widest possible range of bond lengths, the state-averaged CASSCF procedure was employed with density matrix averaging over 22 doublet (six Σ + , seven Π, five ∆, two Φ, and two Σ − ) and 9 quartet (two Σ + , three Π, two ∆, one Φ, and one Σ − ) states, with equal weights for each of the roots. The active space included 7 electrons distributed in 13 orbitals (6 a 1 , 3 b 1 , 3b 2 , 1 a 2 ) that had predominantly O 2p and Y 4d, 5s, 5p, and 6s character; all lower energy orbitals were constrained to be doubly occupied. All valence electrons (4d, 5s Y; 2s, 2p O) were included in the MRCI correlation treatment.
Potential energy, spin-orbit coupling, and dipole moment curves, as well as electronic angular momentum and transition dipole matrix elements were obtained at the MRCI level for the six lowest doublet states. Moreover, the potential curves were calculated using the extended multi-state complete active space secondorder perturbation theory, 73 XMS-CASPT2, with the basis sets augcc-pwCVTZ-PP 71 on Y and aug-cc-pwCVTZ 74 on O (henceforth abbreviated as awCVTZ). In the respective SA-CASSCF calculations, the (7e,13o) active space was employed together with averaging over the lowest six doublet states. In order to remedy issues pertaining to intruder states, a level shift of 0.4 and an IPEA (ionisation potential, electron affinity) shift of 0.5 were employed for XMS-CASPT2.
To calculate the molecular Ω states and respective spin-orbit curves, we used the spin-orbit -MRCI state-interacting approach 75 : the spin-coupled eigenstates were obtained by diagonalizing H es + H SO in a basis of MRCI eigenstates of electrostatic Hamiltonian H es . The matrix elements of H SO were constructed using the one-electron spin-orbit operator accompanying the yttrium pseudopotential.
Spin-orbit effects were also treated more rigorously in relativistic four-component (4c) all-electron calculations employing a Gaussian nuclear model and an accurate approximation to the full Dirac-Coulomb Hamiltonian. 76 The respective spin-free results were obtained with the spin-free Hamiltonian of Dyall. 77 In these calculations, the relativistic TZ-quality basis sets of Dyall 78,79 were used for the Y and O atoms (hereafter referred to as TZ D ). The basis sets were kept uncontracted to provide sufficient flexibility. Electron correlation was taken into account via the equationof-motion CCSD (EOM-CCSD) method 80 with the Y outer-core (4s and 4p) electrons correlated together with the valence electrons. The EOM-EA scheme (adding 1 electron to the closed shell) was applied with the reference defined by the YO + cation and the active space comprising 12 spinors (Y 5s and 4d). For the YO electronic states inaccessible via the EOM-EA procedure, we employed the EOM-IP scheme (removing 1 electron from the closed shell) with the YO − (Y 5s 2 ) anion taken as the reference and an active space composed of 8 spinors (Y 5s and O 2p). The virtual orbital space was truncated by deleting all virtual spinors with orbital energies larger than 15 a.u. In the relativistic calculations of dipole moments, a finite-field perturbation scheme was employed by adding the z-dipole moment operator as a small perturbation to the Hamiltonian. Perturbations with electric field strengths of ±0.0005 a.u. were applied.
The atomic spin-orbit corrections, ∆E SO , utilized in the calculations of the YO atomization energy were obtained from the experimental J-averaged zero-field splittings of the ground state atomic terms 81 : ∆E SO = −77.975 cm -1 (O) and −318.216 cm -1 (Y).
The most sophisticated PEC computations were performed at the coupled cluster singles, doubles, and perturbative triples, CCSD(T), level of theory 82 with a restricted open-shell Hartree-Fock reference and with an allowance for a small amount of spin contamination in the solution of the CCSD equations, i.e., RHF-UCCSD(T). 83,84 Symmetry equivalencing of the ROHF orbitals was performed for the degenerate atomic and molecular electronic states. Both valence (4d, 5s Y; 2s, 2p O) and outer-core (4s, 4p Y; 1s O) electrons were correlated. Scalar relativistic effects were treated with the yttrium pseudopotential described above. Sequences of aug-cc-pwCVnZ-PP 71 (n = T, Q, 5) basis sets for Y were used in conjunction with the corresponding all-electron basis sets aug-cc-pwCVnZ 74 for the O atom. These combinations of basis sets are denoted below as awCVTZ, awCVQZ, and awCV5Z, respectively.
For each point in a grid of r(Y-O) bond lengths, the CCSD(T) calculated energies were extrapolated to the complete basis set (CBS) limit. Three extrapolation schemes were employed. First, a two-point extrapolation of total energies was performed using the formula 85 : where n = 4 and 5 for the awCVQZ and awCV5Z basis sets. This scheme is denoted as CBS1. Second, we employed alternative twopoint (Q5) extrapolations of the Hartree-Fock and correlation energy components. These implied using Eq. 1 for the correlation part and the Karton and Martin formula 86 for the HF energy: this is denoted as CBS2. Third, the CBS estimates were also obtained using the awCVTZ, awCVQZ, and awCV5Z total energies via the three-parameter, mixed Gaussian/exponential expression 87 : where n = 3, 4 and 5 for the awCVTZ, awCVQZ and awCV5Z basis sets, respectively. This is denoted as CBS3. The spectroscopic constants r e , ω e , ω e x e , and α e of YO were obtained from a conventional Dunham analysis 88 using polynomial fits of total energies for bond lengths in the vicinity of the minimum for a given electronic state.
The CCSD(T) calculations of the equilibrium dipole moments, µ e , for a few lowest states were carried out at the corresponding CCSD(T)/CBS1 equilibrium bond lengths. The dipole moments were computed by numerical differentiation of the total energy in the presence of a weak electric field. Finite perturbations with electric field strengths of ±0.0025 a.u. were applied. Since hierarchical sequences of basis sets have been used, the dipole moments were also extrapolated to the CBS limit using the three extrapolation schemes described above.
Most of the ab initio calculations were carried out using the MOLPRO electronic structure package. 89 The relativistic 4c-EOM-CCSD calculations were done with the use of the DIRAC program. 90

Duo calculations
We use the program DUO 65,91 to solve the coupled Schrödinger equation for 6 lowest electronic states of YO. DUO is a variational program capable of solving rovibronic problems for a general (open-shell) diatomic molecule with an arbitrary number of couplings, see, for example, Refs. 58,92-94. All ab initio couplings between these 6 states are taken into account as described below.
The goal of this paper is to provide a qualitative simulation of the electronic spectra of YO based on the ab initio curves. We therefore do not attempt a systematic refinement of the ab initio curves by fitting to the experiment, which will be the subject of future work. In order to facilitate the comparison with the experimental data, we, however, perform some shifts of the T e values and simple scaling of the SOCs (see below).
In DUO calculations, the coupled Schrödinger equation was solved on an equidistant grid of 301 bond lengths r i ranging from r = 1.2 to 3 Å using the sinc DVR method. Our ab initio curves are represented by sparser and less extended grids (see below). For the bond length values r i overlapping with the ab initio ranges, the ab initio curves were projected onto the denser DUO grid using the cubic spline interpolation. The PECs outside the ab initio range were reconstructed using the standard Morse potential form For other curves the following function forms were used 65 : for the short range and for the long range, where A and B are stitching parameters. The vibrational basis set was taken as eigensolutions of the six uncoupled 1D problems for each PEC. The corresponding basis set constructed from 6×301 eigenfunctions was then contracted to include 60 lowest (in terms of the vibrational energy) X functions and 20 from each other state (160 in total). These vibrational basis functions were then combined with the spherical harmonics for the rotational and electronic spin basis set functions. All calculations were performed for 89 Y 16 O using atomic masses.
This study is the first where a DUO calculation has been performed for a system with avoided crossings between curves of the same symmetry. The current version of DUO does not allow for non-adiabatic couplings, and therefore these were ignored in this study. However, despite the expectation that the non-adiabatic coupling effects should be relatively important in the regions near an avoided curve crossing, as we show below, our model neglecting these effects is justified by close agreement with the available experimental spectra.

Electronic structure and potential energy curves of the YO molecule
An overview of the CASSCF PECs for all doublet and quartet states included in the SA-CASSCF procedure is provided in Fig. 1. In the vicinity of the ground state minimum (at ∼1.8 Å), the lowest six CASSCF states are doublets, whereas the quartet states lie at ∼30000 cm -1 above the ground state. The lowest six doublet MRCI PECs ( Fig. 2. For the states X 2 Σ + , A 2 ∆, A 2 Π, and B 2 Σ + , the PECs were obtained in the full bond length range amenable to the underlying CASSCF procedure, 1.58 ≤ r ≤ 2.36 Å, while the C 2 Π and D 2 Σ + curves were calculated through r = 2.04 Å and 1.93 Å, respectively. Extending the MRCI curves for the two upper states beyond those distances would require requesting a greater number of states (while exceeding 4 in a single irreducible representation with the chosen active space would make the MRCI computation unfeasible on the hardware used) or selecting the order of the states in the initial internal CI (e.g., 1, 2, 3, 5 rather than 1, 2, 3, 4), leaving some of them out in each MRCI point. This appeared to affect the smoothness of the resulting PECs. Therefore, we refrained from further pursuing the computations with the same number of MRCI roots in the entire bond length range and simply reduced the number of requested states for longer internuclear distances. For all six doublet states, the XMS-CASPT2 PECs could be obtained in the range 1.59 ≤ r ≤ 2.165 Å, Fig. 3; at larger bond lengths the underlying CASSCF procedure failed to converge. As can be seen from Figs. 2 and 3, the A 2 Π and C 2 Π curves feature an avoided crossing at bond lengths around 2 Å, as do the B 2 Σ + and D 2 Σ + curves in approximately the same region, see Fig. 3. The avoided crossings of both pairs of curves are also seen in the EOM-IP-CCSD calculated PECs, Fig. 4, albeit at shorter distances (1.8 -1.9 Å).
In order to provide deeper insight into the electronic structure of YO, we have analyzed the dominant configurations (configuration state functions) in the MRCI wave functions for the lowest electronic states, Table 1, and the leading atomic orbital contributions in the respective molecular orbitals, Table 2. The analysis indicates that the X 2 Σ + ground state consists mainly of ...10σ 2 11σ 2 5π 4 12σ 1 electron configuration. The prin-cipal contribution to the singly occupied 12σ MO comes from the yttrium 5s atomic orbital, with an admixture of the 5p AO particularly noticeable at longer internuclear distances. The three lowest active MOs, 10σ , 11σ , and 5π +(−) , primarily consist of the oxygen 2s and 2p orbitals whose contributions increase with the bond stretching. Therefore, the bonding in the X 2 Σ + YO ground state can be roughly described as ionic, Y 2+ O 2− , however, with significant covalent character mainly owing to an appreciable participation of the yttrium 4p σ AO in the 10σ MO. Upon the Y-O bond stretching, there is a rapid decrease in the 4p contribution (see Table 2) reducing the covalency and resulting in a concomitant increase in the magnitude of the electric dipole moment in the YO ground state (see below). The first excited state, A 2 ∆, mainly consists of ...10σ 2 11σ 2 5π 4 2δ 1 configuration with the 2δ MO clearly assigned to the Y 4d δ orbital. At shorter bond lengths this state can be reasonably described by the single electron excitation from the ground state, 5s 1 → 4d δ 1 . Upon bond elongation, the weight of the main configuration gradually decreases, approaching ∼50% at bond lengths of about 2.2 Å, whereas at r > ∼2 Å, a few other large-weight configurations emerge, e.g., the three-open-shell ...10σ 2 11σ α 12σ β 5π 4 2δ α configuration (with a weight of 22% at 2.19 Å) that corresponds to the Y + O − bonding (Y 5s 1 4d 1 , O 2p 5 ).

Figure 4
The avoided crossing regions of the EOM-IP-CCSD/TZ D spinfree potential energy curves for the A 2 Π, B 2 Σ + , C 2 Π, and D 2 Σ + electronic states of YO.
As the Y-O distance approaches the avoided crossing point from below, the principal configurations of the A 2 Π and B 2 Σ + states change to those specified above for the C 2 Π and D 2 Σ + states, respectively, while, vice versa, the principal configurations of the C 2 Π and D 2 Σ + states turn into those being the main ones for the A 2 Π and B 2 Σ + states at shorter bond lengths. This alteration of main configurations describes an oxygen-to-metal charge back- Specifically, the avoided crossing point, r ac , chosen to be the point of closest approach of two curves, for the A 2 Π and C 2 Π states amounts to 2.046 Å (XMS-CASPT2) and 1.994 Å (MRCI), with the energy gap, ∆E ac , of 366 cm -1 and 243 cm -1 , respectively. For the XMS-CASPT2 B 2 Σ + and D 2 Σ + curves, r ac = 2.064 Å and ∆E ac = 1484 cm -1 . Notably, the principal configuration interchange between the B 2 Σ + and D 2 Σ + states occurs at a slightly shorter internuclear distance: ∼1.95 Å (XMS-CASPT2), ∼1.92 Å (MRCI). At the EOM-IP-CCSD level, the avoided crossing characteristics are: r ac = 1.911 Å, ∆E ac = 231 cm -1 for A 2 Π and C 2 Π curves, and r ac = 1.832 Å, ∆E ac = 272 cm -1 for B 2 Σ + and D 2 Σ + curves.
The CCSD(T) calculations were carried out for the six lowest doublet states. The reference configurations for each state were as follows: The CCSD(T) energies were obtained in the ranges: At longer distances (as well as shorter ones for C 2 Π and D 2 Σ + ), the coupled-cluster calculations failed due to severe CCSD convergence problems. For the X 2 Σ + , A 2 ∆, A 2 Π, and B 2 Σ + states, the distances shorter than 1.0 Å were not considered because relative energies of excited states already exceed 350000 cm -1 at this point. The respective CCSD(T)/CBS1 PECs are shown in Fig. 5. Since single reference methods are not suitable for describing avoided crossings, e.g., those between the PECs of the A 2 Π and C 2 Π, and B 2 Σ + and D 2 Σ + states, the CCSD(T) calculated PECs can be viewed as corresponding to the diabatic presentation of these states. It is clearly seen that the CCSD(T) A 2 Π and C 2 Π curves cross each other at 2.031 Å, as do the B 2 Σ + and D 2 Σ + ones at 2.013 Å.  Table 3 summarizes the optimized bond lengths r e , harmonic vibrational frequencies ω e , equilibrium dipole moments µ e , and adiabatic excitation energies T e of the low-lying doublet states calculated at the XMS-CASPT2, MRCI and EOM-CCSD levels as well as the results from earlier theoretical studies. 37,53 Our data in the column entitled "C 2 Π" were obtained from the second minimum in the adiabatic PEC of the A 2 Π state, i.e., they can be ascribed to the diabatic representation of the A 2 Π and C 2 Π states.

Spectroscopic constants and excitation energies
The results of our EOM-CCSD calculations indicate that the anionic reference is less suitable for describing YO than the cationic one. In Table 3, more accurate cationic-reference EOM-EA-CCSD spectroscopic constants are listed for all states except for the C 2 Π and D 2 Σ + ones which are not accessible via the electron attachment procedure and therefore were described at the EOM-CCSD level only via EOM-IP.
The results given in Table 3 are obviously inferior to those obtained from high-level CCSD(T) calculations including corevalence correlation and extrapolation to the CBS limit. The CCSD(T) results are collected in Table 4 together with the experimental data available to date. The spread in the CCSD(T)/CBS results from different CBS extrapolation schemes serve as a rough estimate of the uncertainty in extrapolation. The good agreement between the CBS estimates and experimentally determined spectroscopic properties of the X 2 Σ + , A 2 ∆, A 2 Π, and B 2 Σ + electronic states demonstrates the high accuracy in the CCSD(T)/CBS PECs of these states for bond lengths in the vicinity of the PECs minima, and is indicative of the mild MR character of the respective electronic wave functions.

Figure 6
The CCSD/awCV5Z calculated T 1 diagnostics of YO An insight into the reliability of the CCSD(T) PECs over the entire bond length range explored, and for all electronic states considered, including those not yet characterized experimentally, can be provided by using the MR diagnostic criteria commonly employed to identify the suitability of single reference wavefunctionbased methods: T 1 , 95 the Frobenius norm of the coupled cluster amplitude vector related to single excitations, and D 1 , 96 the matrix norm of the coupled cluster amplitude vector arising from coupled cluster calculations. The utility of different MR diagnostics has been examined in earlier studies 97,98 on various 3d and 4d TM species. The following criteria have been proposed 98 as a gauge for the latter to predict the possible need to employ multireference wavefunction-based methods while describing energetic and spectroscopic molecular properties:  The CCSD/awCV5Z T 1 plots vs. Y-O bond length are shown in Fig. 6. The similar D 1 plots are illustrated in Fig. S1 of the supplementary material. At shorter bond lengths, the diagnostics amount to 0.02-0.03 (T 1 ) and 0.05-0.12 (D 1 ), remaining below the MR thresholds down to 1.4 Å for most states. Upon bond stretch, T 1 and D 1 rapidly increase, typically exceeding the MR threshold at 2.1-2.2 Å. The behaviour of these diagnostics for the C 2 Π state is a Table 4 CCSD(T) spin-free equilibrium constants of YO in its low-lying electronic states: The dissociation energy D 0 (eV) referring to the ground electronic state X 2 Σ + , the excitation energies T e (cm -1 ) of the A 2 ∆, A 2 Π, B 2 Σ + , C 2 Π D 2 Σ + , and a 4 Π states, bond length r e (Å), spectroscopic constants ω e (cm -1 ), ω e x e (cm -1 ) and α e (cm -1 ), and dipole moment µ e (D). notable exception: the numerical values of both T 1 and D 1 remain well below the MR threshold throughout the bond length range studied. The D 2 Σ + state is also noteworthy: its T 1 and D 1 diagnostics are indicative of the CCSD D 2 Σ + wave function retaining its SR character in much narrower range of bond lengths compared to the other doublet states under study.
The relative importance of SR/MR character of YO can also be guessed with the use of spin contamination appearing in RHF-UCCSD calculations as a result of unrestricted spin at the CCSD level. According to Jiang et al., 97 spin contamination with <S 2 − S 2 z − S z > greater than 0.1 in an RHF-UCCSD wave function can be viewed as a strong indication of nondynamical correlation in an open-shell system. Plotting spin contamination vs. bond length, Fig. 7, clearly indicates the severe admixture of higher spin states in the CCSD A 2 ∆, A 2 Π, B 2 Σ + and D 2 Σ + wave functions at bond lengths beyond 2.2-2.3 Å. Greater extent of spin contamination at longer internuclear distances can obviously be associated with larger values of T 1 and D 1 (Fig. 6 and Fig. S1 of the supplementary material) exceeding the established MR thresholds. It is instructive to compare the MR diagnostics discussed above with the weights of the principal configurations, C 2 0 , in the MRCI wavefunctions of YO (see Table 1). At shorter bond lengths, the C 2 0 values amount to ∼0.73 for the B 2 Σ + state and 0.79 -0.83 for the remaining doublet electronic states under study. These values are smaller than the threshold, C 2 0 ≥ 0.90, proposed in Refs. 97,98 to recognize the wave function strongly dominated by a single configuration. It should, however, be noted that this criterion was established by analyzing the CASSCF wavefunctions, whereas the C 2 0 of the entire MRCI wavefunction differs from that of the CASSCF reference function due to the contributions of external configurations, which make C 2 0 a smaller number. Upon the YO bond stretch, there is a gradual decrease in the weights of the configurations serving as a reference for the coupled cluster treatment of the X 2 Σ + , A 2 ∆, and A 2 Π states. This indicates greater multireference character of the respective wave functions at longer bond distances, as do the CC-based MR diagnostics. The reference configuration for the C 2 Π state has approximately the same weight, C 2 0 ∼ = 0.83, in the MRCI wavefunctions of YO throughout the bond length range explored, behaving just like the respective CC-based MR diagnostics. These examples of the CCSD -MRCI correlations imply that the CC-based MR diagnostics can be capable of providing qualitative data about the relative accuracy in the single-reference coupled cluster calculation results not only for near-equilibrium regions of electronic ground states, but also for excited states in a wider range of molecular geometries.
In general, the present analysis indicates essentially single reference character of the YO low-lying doublet states over most part of bond length range explored in our work and hence high accuracy in the respective domains of the CCSD(T) PECs. It may also be indicative of accuracy degradation at larger bond lengths, implying the need for additional adjustments of the CCSD(T) PECs. Nevertheless, the bond length range associated with high-energy sections of PECs is expected to have a limited impact on the simulated spectra.

Quartet states
We have studied five low-lying quartet electronic states of YO at the CASSCF, CASPT2, CASPT3, and MRCI levels of theory using the aVTZ basis set. The results are shown in Table 5 together with the earlier theoretical findings. 53 The lowest quartet, a 4 Π, was also studied at the CCSD(T) level (Table 4). At larger internuclear distances, e.g., at r = 2.19 Å, all the quartets feature similar orbital Table 5 Theoretical spin-free equilibrium constants of YO in its low-lying quartet states: adiabatic excitation energy T e -cm -1 , bond length r e -Å, vibrational frequency ω e -cm -1 . The aVTZ basis set has been used throughout. character corresponding to the Y 5s 1 4d 1 , O 2p 5 electron configuration consistent with the Y + O − bonding (see Table 6).
The results for the YO quartet states obtained in our work agree with those of Langhoff and Bauschlicher 53 , Table 5, except for the symmetry of the lowest quartet state that was reported 53 to be 4 Φ rather than 4 Π.
The single-reference CCSD(T) method is expected to yield quite accurate results for the a 4 Π state of YO since its MR diagnostics, C 2 0 = 0.90, T 1 = 0.024, and D 1 = 0.078, indicate essentially SR character of the a 4 Π wave function in the vicinity of the a 4 Π PEC minimum, 2.00 -2.25 Å. The very large CCSD(T)/CBS excitation energy of the a 4 Π state, 29600 cm -1 , suggests that the quartet states in YO are too high in energy to significantly affect the spectroscopy of its low-lying doublet states.

SO coupling
Spin-orbit coupling effects were studied in a perturbative fashion at the MRCI level and more rigorously at the 4c-EOM-CCSD level of theory including spin from the outset. The 4c-EOM-CCSD calculated spin-orbit coupling effects on the spectroscopic constants of YO are shown in Table 7. The theoretical spin-orbit coupling constants, SOCCs (A SO ), can be compared with the relevant experimental numbers for the A 2 ∆ and A 2 Π electronic states of YO reported previously. 11,13,38 The A SO (A 2 ∆) SOCCs of 336 cm -1 and 313 cm -1 obtained at the MRCI SI-SO and 4c-EOM-CCSD levels, respectively, agree well with each other and with the experimental number of 339 cm -1 determined by Chalek and Gole 38 . The calculation results for A SO (A 2 Π), 346 cm -1 (MRCI SI-SO) and 438 cm -1 (4c-EOM-CCSD), are also in reasonable agreement with the respective experimental value of 428 cm -1 . 11,13 As the Y-O distance reaches the avoided crossing point between A 2 Π and C 2 Π, Table 6 Main configurations in the low-lying quartet electronic states of YO derived from analyzing the MRCI/aVTZ wave function at a bond length of 2.19 Å. Weight of each configuration is ∼45%. Table 1 for designations. the A SO (C 2 Π) and A SO (A 2 Π) values change their sign: the A 2 Π 3/2 spin component of the A 2 Π state becomes lower in energy than its A 2 Π 1/2 counterpart, and vice versa for the spin-coupled components of the C 2 Π state (see Fig. 8). There is also a change in the absolute values of A SO (A 2 Π) and A SO (C 2 Π): at r < r ac , |A SO (A 2 Π)| is much lower in magnitude than |A SO (C 2 Π)| and vice versa at r > r ac . However, the numerical values of A SO (A 2 Π) at r > r ac determined with the MRCI SI-SO and 4c-EOM-CCSD methods, e.g., −45 cm -1 and −186 cm -1 , respectively, at a bond length of 2.04 Å, are in less satisfactory agreement with each other than those at r < r ac .
The SOC matrix elements between various doublet states of YO, which also accurately account for the corresponding phases, are shown in Fig. 9 as a function of r(Y-O). The relative phases of the couplings are important when used for solving the nuclear motion problem as part of the coupled Schrödinger equation, see, for example, discussion by Patrascu et al. 92 The full details of the ab initio coupling curves including the magnetic quantum numbers are provided as part of the supplementary data.

Dipole moment, transition dipole moment, and electronic angular momentum curves of YO
The CCSD(T)/CBS dipole moment of 4.61 D for the X 2 Σ + state of YO (  45 from the more precise microwave measurement, 4.524 (7) 27 , the reason for µ(A 2 Π 3/2 ) being larger than µ(A 2 Π 1/2 ) in YO in contrast to ScO is the smaller energy gap of the A 2 Π and A 2 ∆ states, which results in mixing between the Ω = 3/2 spin-orbit components of these states. This idea is, however, based on low-level ab initio computations by Langhoff and Bauschlicher 53 , which predicted the A 2 ∆ state in YO to lie 200 cm -1 higher than A 2 Π, whereas for ScO the analogous calculations 102 yielded a difference of about 1900 cm -1 , with A 2 ∆ being lower in energy. In fact, the A 2 ∆ state lies around 1800 cm -1 lower than A 2 Π in YO and 1500 cm -1 lower in ScO, as evidenced by experimental works of Chalek and Gole 38,49 , i.e., the A 2 Π -A 2 ∆ energy gap in YO exceeds that in ScO. Also, at high levels of theory including SO coupling, the PECs for the Ω = 3/2 components of A 2 Π and A 2 ∆ lie quite far apart (see the excitation energies in Table 7), and their mixing is almost negligible. Furthermore, it is worth noting that for another analogous molecule, LaO, the experimental data 103 also indicate that µ e (A 2 Π 1/2 ) > µ e (A 2 Π 3/2 ). To shed more light on the alleged different order of the dipole moment values for the spin-orbit components of the A 2 Π state in YO compared to ScO and LaO, we performed additional 4c-EOM- Since the 4c-EOM-EA-CCSD dipole moments are expected to be overestimated by at least 0.5 D, one can consider the theoretical results to be in reasonable agreement with experiment. In light of our results, the experimental dipole moments 27 for the two Ω components of the A 2 Π state of YO need to be revisited.
The MRCI DMCs and TDMCs of YO are shown in Fig. 10. The EAMCs are shown in Fig. S2 of the supplementary material. All these curves as well as the SOC ones (Fig. 9) exhibit irregular behaviour at bond lengths around r ∼2 Å due to strong changes in the A 2 Π, B 2 Σ + , C 2 Π and D 2 Σ + wave functions over the avoided crossing region.

Results of DUO calculations
For DUO calculations, we selected the following set of curves representing our highest level of theory: the CCSD(T)/CBS PECs shown in Fig. 5 and MRCI SOCs, (T)DMCs, and EAMCs shown in Figs. 9 and 10, as well as Fig. S2 of the supplementary material. Due to limitations of single reference CCSD(T), the CCSD(T) curves for the A 2 Π, C 2 Π and B 2 Σ + , D 2 Σ + states do not exhibit avoided crossings and hence are not consistent with the MRCI property curves. To alleviate this deficiency, these four PECs were transformed by simply switching the corresponding energy points between A and C ( 2 Π states) as well as those between B and D ( 2 Σ + states) at r > r ac . We have decided to apply this rather simplistic procedure because it has marginal effect on the overall accuracy of our model and is sufficient for the goal of this pure ab initio study, not aiming at spectroscopic accuracy. A proper diabatic representation of the YO electronic states will be, however, important when refining the ab initio curves by fitting to experiment, 93 which is a goal of future work. In this study we work directly with the ab initio data in the grid representation without representing ab initio curves analytically. We do not perform any diabatizations here, which is often useful for representing the variation of the data with respect to the bond length in an intuitive and more compact form. One of the artifacts of this choice to use the ab initio curves directly is that the crossing points of the MRCI PECs hence the points of drastic change in the MRCI property curves differ by a few hundredths of Å from the crossing points of the CCSD(T)/CBS PECs. Again, this has small impact on the overall agreement of the current model with the experiment. However, a more accurate study will require a more consistent treatment of the crossing points. Our preferred choice would be to use the CCSD(T)/CBS values of the corresponding crossing points. Table 7 4c-EOM-CCSD/TZ D molecular properties of YO in its low-lying spin-coupled electronic states (EOM-IP for C 2 Π and EOM-EA for the remaining states): bond length r e -Å, vibrational frequency ω e , adiabatic excitation energy T e -cm -1 . The respective spin-orbit effects, ∆ SO , are provided as well. −0.0 21994 +98 Figure 10 MRCI/aVTZ electric dipole moment curves of YO: diagonal (above), non-diagonal µ x (middle) and non-diagonal µ z (below).
The DUO rovibronic wavefunctions of YO in conjunction with the ab initio TDMCs were then used to produce Einstein A coef-ficients for all rovibronic transitions between states considered in this work covering the wavenumber range from 0 to 40000 cm -1 and J ≤ 180.5. These Einstein A coefficients and the corresponding energies from the lower and upper states involved in each transition were organized as a line list using the ExoMol format. 106 This format uses a two file structure with the energies included into the States file (.states) and Einstein coefficients appearing in the Transitions file (.trans). This ab initio line list is available from www.exomol.com. The ExoMol format has the advantage of being compact and compatible with our intensity simulation program EXOCROSS 107 (see below).

Vibronic energies
In Table 8 we compare our computed vibrational excitations at J = 0.5 and J = 1.5 (as proxy for vibrational band centres) of 89 Y 16 O with the experimentally derived values. Based on this comparison, as an ad hoc improvement we applied the following shifts to PECs of the excited states: +9.509 cm -1 (A 2 Π), +81.096 cm -1 (A 2 ∆), −134.301 cm -1 (B 2 Σ + ), and +358.626 cm -1 (D 2 Σ + ). These shifts are small compared to the observed minus calculated differences often encountered in calculations of electronic term values for transition metal oxides. 60 We also scaled the SOC of A 2 Π by 1.14 in order to increase the SO splitting of v = 0 by about 33 cm -1 . Even though we are not targeting fully quantitative accuracy in this work, without such empirical shifts it would be difficult to reproduce band heads even qualitatively.
To allow for a direct comparison with the observed spectra of 89 Y 16 O, we generated a line list covering rotational excitations up to J = 190 and the energy/wavenumber range up to 40,000 cm -1 , with a lower state energy cutoff of 16,000 cm -1 .

Partition function
The partition function of 89 Y 16 O computed using our ab initio line list is shown in Fig. 11, which is compared to that recently reported by Barklem and Collet 110 . Since 89 Y has a nuclear spin degeneracy of two, we have multiplied Barklem and Collet's partition function by a factor of two to compensate for the different conventions used; we follow HITRAN 111 and include the full nuclear spin in our partition functions. The partition function of YO was also reported by Vardya 112 , which is shown in Fig. 11. All three partition functions are almost identical for their ranges of validity.

Spectral comparisons
Using the ab initio 89 Y 16 O line list, spectral simulations were performed with our code EXOCROSS. 107 EXOCROSS is an open source Fortran 2003 code with the primary use to produce spectra of molecules at different temperatures and pressures in the form of cross sections using molecular line lists as input.  Here we use the YO line list generated with DUO in the Ex-oMol format, the description of which can be found, e.g., in Yurchenko et al. 107 or Tennyson et al. 106 . EXOCROSS can be accessed via http://exomol.com/software/ or directly at https://github.com/exomol. Amongst other features, EX-OCROSS can generate spectra for non-local thermal equilibrium conditions characterized with different vibrational and rotational temperatures, lifetimes, Landé g-factors, partition and cooling functions.
An overview of the YO absorption spectra in the form of cross sections at the temperature T = 2000 K is illustrated in Fig. 12.
Here, a Gaussian line profile with a half-width-at-half-maximum (HWHM) of 5 cm -1 was used. This figure shows contributions from each electronic band originating from the ground electronic state. The strongest bands are A 2 Π-X 2 Σ + and B 2 Σ + -X 2 Σ + . The visible A-X band is known to be important for the spectroscopy of cool stars. The C state is of the same symmetry as A, however, the corresponding band C-X is much weaker due to the small Franck-Condon effects. The A 2 ∆-X 2 Σ + band is forbidden and barely seen in Fig. 12, however, it is strong enough to be experimentally known. 39 Figure 13 shows a simulated emission spectrum of the strongest orange system YO (A 2 Π-X 2 Σ + , (0,0)), which is compared to the experiment of Badie and Granier 31 (from the plume emission close to the liquid Y 2 O 3 surface). It is remarkable that even pure ab initio calculations (after modest adjustment of the corresponding T e value by +9.509 cm -1 ) provide very close reproduction of experiment. It shows that our line list at the current, ab initio quality should be useful for modelling spectroscopy of exoplanets and cool stars in the visible region.  Figure 14 illustrates the A 2 ∆ -X 2 Σ + (0,0) forbidden band in emission simulated for T = 77 K compared to the experimental spectrum of Simard et al. 39 . Here, a shift of +81.096 cm -1 was applied to the T e value of the A 2 ∆ state. In spectral simulations, this region appears to be contaminated by the dipole-allowed hot A-X transitions, which are not necessarily very accurate in this region. We therefore applied a filter to select the A 2 ∆ -X 2 Σ + transitions only. The difference in shape of the spectra can be attributed either to the non-LTE (Local Thermal Equilibrium) effects present in the experiment or broadening effects, which we have not attempted to model properly. This figure is only to illustrate the generally good agreement of the positions of the rovibronic lines in this band.  Figure 15 shows a series of absorption bands compared to the measurements of Zhang et al. 37 who observed bands in both the B 2 Σ + -X 2 Σ + and D 2 Σ + -X 2 Σ + systems in a heavily nonthermal environment where the vibrations were hot and the rotations cooled to liquid nitrogen temperatures. In this case of multiband system it was important to include at least some non-LTE effects by treating it using two temperatures, vibrational and rotational, assuming that the corresponding degrees of freedom are in LTE. The rotational temperature T rot = 77 K was set to value specified by Zhang et al. 37 , while the vibrational temperature was adjusted to T vib = 2000 K to better reproduce the experimental spectrum. The spectrum is divided into five spectroscopic windows (I-V) which are also detailed in Table 9. In order to match the positions of the vibronic bands in the experiment, some of the windows were shifted. For example, the D 2 Σ + -X 2 Σ + (1,0) band was shifted by about 76.5 cm -1 . This shift is an indication of the inaccuracy of our model to reproduce the vibrationally excited states of D 2 Σ + . This is not surprising considering the complexity of the quantum-chemistry part of these systems as well as of the nuclear motion part. The avoided crossing with the B 2 Σ + state leads to very complex shapes of the D 2 Σ + PEC and of the SO and electronic angular momentum coupling curves with the A and C states. The corresponding SOCs of the B and D states with the nearby state C are also relatively large, ∼ 30 cm -1 and 80 cm -1 , respectively (see Fig. 9), and therefore important. Besides, the D PEC is rather shallow with the equilibrium in the vicinity of the avoided crossing point, which also complicates the solution. An accurate description of the B and D curves would require diabatic representations before attempting any empirical refinement by fitting to the experiment. In all cases our simulations, while not perfect, show striking agreement with the observed spectra. Table 9 Five spectroscopic windows (cm -1 ) used to compare five vibronic bands of YO (B and D) in Fig. 15. Experiment is by Zhang et al. 37  Comparison of our computed emission spectra to the measurements of Zhang et al. 37 Our simulations assumed a cold rotational temperature of T rot = 77 K and a hot vibrational temperature of T vib = 2000 K. The Gaussian line profile of HWHM = 0.1 cm -1 was used.

Lifetimes
The lifetimes of 89 Y 16 O in the A 2 Π and B 2 Σ + states (v ≤ 2) were measured by Liu and Parson 12 using laser fluorescence detection of nascent product state distributions in the reactions of Y with O 2 , NO, and SO 2 . Some lifetimes were also measured by Zhang et al. 37 and computed ab initio by Langhoff and Bauschlicher 53 . Table 10 presents a comparison of these results with our calculations with DUO, 50 showing value of the states with corresponding lowest J and the positive parity. It can be seen that our A 2 Π state lifetimes appear to be shorter than the observed ones. This suggests that the A 2 Π -X 2 Σ + transition dipoles may be slightly too large. Good agreement is obtained for the lifetimes of the B 2 Σ + states, while the D 2 Σ + state lifetimes are underestimated by a factor of 2 indicating that the corresponding transition dipole moments D-B and D-X, or at least one of them might be too large.

Conclusion
In this work, a composite approach to accurate first-principles description of the spectroscopy of open-shell TM-containing diatomics is proposed and its high efficiency is demonstrated taking the example of the yttrium oxide molecule. The approach is based on the combined use of single reference coupled cluster and multireference methods of electronic structure theory, accompanied with a thorough joint analysis of the SR/MR character of the molecular wave function. A full set of potential energy, (transition) dipole moment, spin-orbit, and electronic angular momenta curves for the lowest 6 electronic states of YO was produced ab initio using a combination of the CCSD(T)/CBS and MRCI methods. These curves were then used to solve the fully coupled Schrödinger equation for the nuclear motion using the DUO program. Given the complexity of the system under study, the results show remarkably good agreement with the experiment. Our ultimate goal is to produce an accurate, empirical line list for 89 Y 16 O for applications in modelling the spectroscopy of atmospheres of exoplanets and cool stars. This will require a refinement of the ab initio curves by fitting to the experimental data in the diabatic representation as well as inclusion of the non-adiabatic coupling effects and will be addressed in future work. The A 2 Π band of YO has strong absorption in the visible region, i.e. where the stellar radiation usually peaks. Such systems are known to cause the temperature inversion in atmospheres of exoplanets, similar to the inversion caused by TiO and VO in giant exoplanets 113 . Opacities of such species are crucial in modelling the degree of temperature inversion in giant exoplanets. YO is yet to be detected in exoplanetary atmospheres and this work is meant to provide the necessary spectroscopic data. YO is one of the few molecules with the strong potential for laser-cooling applications, 18 which have widely ranging applications, from quantum information and chemistry to searches for new fundamental physics. The results of this work will help to model the cooling properties of YO and thus will be important for designing and implementing laser-cooling experiments.
The ab initio curves of YO obtained in this study 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 computed line list can be obtained from www.exomol.com. This is given in the ExoMol format 106 which also includes state-dependent lifetimes. The line list can be directly used with the EXOCROSS program to simulate the spectral properties of YO. UCL Legion High Performance Computing Facility (Legion@UCL), and associated support services, in the completion of this work, along with the Cambridge COSMOS SMP system, part of the STFC DiRAC HPC Facility supported by BIS National E-infrastructure capital grant ST/P002307/1 and ST/R002452/1 and STFC operations grant ST/R00689X/1. DiRAC is part of the National e-Infrastructure. A. N. S. and V. G. S. acknowledge support from the Ministry of Science and Higher Education of the Russian Federation (Project No. 4.3232.2017/4.6).