Diagonal Born-Oppenheimer Corrections to the Ground Electronic State Potential Energy Surfaces of Ozone: Improvement of Ab Initio Vibrational Band Centers for the 16O3, 17O3 and 18O3 Isotopologues

Mass-dependent diagonal Born-Oppenheimer corrections (DBOC) to the ab initio electronic ground state potential energy surface for tseveral isotopologues of the ozone molecule are reported for the first time. The comparison with experimental band centers shows a significant improvement of the accuracy with respect to the best Born-Oppenheimer (BO) ab initio calculations reducing the total root-mean-squares (calculated observed) deviations by about factor of two. For the set of 16O3 vibrations up to five bending and four stretching quanta, the mean (calculated - observed) deviations drop down from 0.7 cm-1 (BO) to about 0.1 cm-1, with the most pronounced improvement seen for bending states and for mixed bend-stretch polyads. In case of bending band centers directly observed under high spectral resolutions, the errors are reduced by more than order of magnitude from observed levels, approaching nearly experimental accuracy. New sets of ab initio vibrational states can be used for improving spectroscopic effective models for analyses of observed high-resolution spectra, particularly in cases of accidental resonances with ,,dark'' states requiring accurate theoretical predictions. Mass-dependent diagonal Born-Oppenheimer corrections (DBOC) to the ab initio electronic ground state potential energy surface for the main 16 O 3 isotopologue and for homogeneous isotopic substitutions 17 O 3 and 18 O 3 of the ozone molecule are reported for the ﬁrst time. The sys-tem being of strongly multiconﬁgurational character, Multireference Conﬁguration Interaction wave function ansätze with different complete active spaces were employed. Reliable DBOC calculations with the targeted accuracy were possible to conduct up to about the half of the dissociation threshold D 0 . The comparison with experimental band centers shows a signiﬁcant improvement of the accuracy with respect to the best Born-Oppenheimer (BO) ab initio calculations reducing the total root-mean-squares (calculated – observed) deviations by about factor of two. For the set of 16 O 3 vibrations up to ﬁve bending and four stretching quanta, the mean (calculated â˘A¸S observed) deviations drop down from 0.7 cm − 1 (BO) to about 0.1 cm − 1 , with the most pronounced improvement seen for bending states and for mixed bend-stretch polyads. In case of bending band centers directly observed under high spectral resolutions, the errors are reduced by more than order of magnitude down to 0.02 cm − 1 from observed levels, approaching nearly experimental accuracy. A similar improvement for heavy isotopologues shows that reported DBOC corrections almost remove the systematic BO errors in vibrational levels below D 0 / 2, though the scatter increases towards higher energies. The possible reasons for this ﬁnding, as well as remaining issues are discussed in detail. The reported results provide an encouraging accuracy validation for the multireference methods of the ab initio theory. New sets of ab initio vibrational states can be

: Comparison of observed and calculated band centers (pure vibration energies) for the 16 O 3 isotopologue of ozone. Each level is given by its symmetry ("SYMM"), vibrational quantum assignment ("v 1 ,v 2 ,v 3 "), observed energy ("OBS"), and a set of calculated energies: Born-Oppenheimer PES without DBOC corrections ("BO"), and DBOC-corrected energy (BO+DBOC(n,n)). All band centers are given in wavenumbers. For each calculation, the zero point energies ("ZPE") are provided.  Table 2: Comparison of observed and calculated band centers (pure vibration energies) for the 17 O 3 isotopologue of ozone. Each level is given by its symmetry ("SYMM"), vibrational quantum assignment ("v 1 ,v 2 ,v 3 "), observed energy ("OBS"), and a set of calculated energies: Born-Oppenheimer PES without DBOC corrections ("BO"), and DBOC-corrected energy (BO+DBOC (6,6)). All band centers are given in wavenumbers. For each calculation, the zero point energies ("ZPE") are provided.  Table 3: Comparison of observed and calculated band centers (pure vibration energies) for the 18 O 3 isotopologue of ozone. Each level is given by its symmetry ("SYMM"), vibrational quantum assignment ("v 1 ,v 2 ,v 3 "), observed energy ("OBS"), and a set of calculated energies: Born-Oppenheimer PES without DBOC corrections ("BO"), and DBOC-corrected energy (BO+DBOC(6,6)). All band centers are given in wavenumbers. For each calculation, the zero point energies ("ZPE") are provided. Mass-dependent diagonal Born-Oppenheimer corrections (DBOC) to the ab initio electronic ground state potential energy surface for the main 16 O 3 isotopologue and for homogeneous isotopic substitutions 17 O 3 and 18 O 3 of the ozone molecule are reported for the first time. The system being of strongly multiconfigurational character, Multireference Configuration Interaction wave function ansätze with different complete active spaces were employed. Reliable DBOC calculations with the targeted accuracy were possible to conduct up to about the half of the dissociation threshold D 0 . The comparison with experimental band centers shows a significant improvement of the accuracy with respect to the best Born-Oppenheimer (BO) ab initio calculations reducing the total root-mean-squares (calculated -observed) deviations by about factor of two. For the set of 16 O 3 vibrations up to five bending and four stretching quanta, the mean (calculated âȂŞ observed) deviations drop down from 0.7 cm −1 (BO) to about 0.1 cm −1 , with the most pronounced improvement seen for bending states and for mixed bend-stretch polyads. In case of bending band centers directly observed under high spectral resolutions, the errors are reduced by more than order of magnitude down to 0.02 cm −1 from observed levels, approaching nearly experimental accuracy. A similar improvement for heavy isotopologues shows that reported DBOC corrections almost remove the systematic BO errors in vibrational levels below D 0 / 2, though the scatter increases towards higher energies. The possible reasons for this finding, as well as remaining issues are discussed in detail. The reported results provide an encouraging accuracy validation for the multireference methods of the ab initio theory. New sets of ab initio vibrational states can be used for improving spectroscopic effective models for analyses of observed high-resolution spectra, particularly in cases of accidental resonances with "dark" states requiring accurate theoretical predictions.

Introduction
Accurate potential energy surfaces (PES) along the nuclear coordinates are pre-requests for reliable theoretical studies of molecular spectroscopy and dynamics. spectra of the ozone molecule is mandatory for spectra analyses because of its importance for atmospheric applications. [1][2][3][4] Since early ab initio works, [5][6][7] it has been recognized that the ozone molecule possesses a quite complicated electronic structure 1,[8][9][10][11] attracting much attention both for the study of the ground [12][13][14] and excited electronic states. 15 On the experimental side, the incentive for in-depth investigations was related to isotopic anomalies in the ozone formation discovered both in the atmosphere and in the laboratory settings. [16][17][18][19][20][21] Over the years this was a motivation for many studies of vibration energy patterns, 12,[22][23][24] infrared spectra of ozone isotopologues (see refs. 25-31, and references therein) and for the dynamics of isotopic exchange reactions . [32][33][34][35][36][37][38] The non-adiabatic coupling of electronic states has been discussed in refs. 39,40. However, to our knowledge, all available full-dimensional PESs of the ozone molecule 1,12-14 have been computed within the framework of the Born-Oppenheimer approximation. The present work aims at first ab initio calculation of the mass-dependent diagonal Born-Oppenheimer corrections (DBOC) for the 3-dimensional ozone PES and their application to vibrational band centers of the 16 O 3 , 17 O 3 and 18 O 3 isotopic species. We show that these contributions permit a significant improvement of vibrational calculations with respect to the most accurate available BO surface.
The paper is structured as follows. The ansatz of DBOC calculations, the wave function model for ozone DBOC and the corresponding computational methodology are described in Section 2. The analytic PES model employed for the fit of DBOC electronic energy corrections in terms of nuclear geometries is discussed in Section 3. The subsequent Sections 4 and 5 are devoted to a detailed study of the mass-dependent corrections to vibrational states of the 16 O 3 , 17 O 3 and 18 O 3 isotopologues and a comparison to experimental data with the conclusions and discussions in Section 6.

Diagonal Born-Oppenheimer Correction
The Diagonal Born-Oppenheimer Correction (DBOC) is the leading correction term to the Born-Oppenheimer Approximation (BOA) first given by Born and Huang: 41 (1) with Ψ denoting the normalized electronic wave function obtained within the BOA, andT N being the nuclear kinetic energy operator. The integration in eqn 1 is done over all electronic coordinates represented by r. The DBOC thus takes into account the dependence of the electronic wave function on the nuclear coordinates R, through the calculation of the nuclear kinetic energy contribution. It has the advantage that the adiabatic picture of the BOA is retained so that the notion of a potential energy surface (PES) is still possible. However, the PES becomes dependent on the nuclear mass and will thus be different for, e.g. various isotopologues of the same molecule.
The ab initio calculation of ∆E DBOC is possible with analytic derivative techniques. Following the pioneer work of Sellers and Pulay 42 and Handy et al. 43 who first presented a formula for the evaluation of ∆E DBOC at the Hartree-Fock Self-Consistent Field (HF-SCF) level, several implementations were reported for correlated electronic wave functions. [44][45][46][47][48][49][50][51][52] The work some of the present authors 51 made first possible to calculate ∆E DBOC at various levels of Configuration Interaction (CI) and Coupled Cluster (CC) theory via analytic techniques. The implementation of their formulae for the single-reference Coupled Cluster Singles and Doubles (CCSD) and Configuration Interaction Singles and Doubles (CISD) models, as well as for the Møllet-Plesset perturbation theory are available today in the CFOUR 53 program. Using an interface to the MRCC 54 program system, calculation of ∆E DBOC is possible at not only any truncation level of singlereference CC and CI theories, but also for various multiconfigurational CI and CC models. [55][56][57][58] The evaluation of ∆E DBOC requires the explicit calculation of first order wave function response pa-rameters with respect to nuclear coordinates, and it is thus computationally quite demanding.
The DBOC is often regarded as a predominantly one-electron effect, with the contribution of electron correlation being small, even negligible. However, as shown in refs. 51 and 52, this is not at all the case if accurate energy differences are desired, to which the account of electron correlation contribution in DBOC falls easily in the range of the total DBOC effect. Treating DBOC at correlated levels is thus necessary in studies aiming at subchemical (kJ/mol or better) accuracy. Fortunately, however, even a relatively low level treatment of dynamic correlation is able to recover a large fraction of the DBOC electron correlation contribution. 52 In terms of the basis set size, is also shown in ref. 51 that the total DBOC can essentially be considered as converged in triple-ζ quality basis sets, but several wavenumbers away from the basis set limit at the double-ζ level. DBOC corrections have been included in accurate PES calculations for several molecules (see refs. 47, 59-64 and references therein), however, to our knowledge have never been accounted for the ozone molecule.

Wave function model for ozone DBOC
In the case of ozone, the complexity of the electronic structure makes the use of a multiconfigurational (MC) wave function model necessary, even in points near the equilibrium structure. 1,9,10,12 The key step in these calculations is the choice of the reference space, which is in most cases a CAS (complete active space) where the "most important" orbitals (the so called active orbitals) are selected and the"most important" electrons (active electrons) are distributed among them. Technically, CAS is denoted by CAS(m,n), where m is the number of active electrons and n is the number of active orbitals. From the chemical point of view, the model that selects all valence orbitals and electrons as active (often referred to as the full-valence CAS space) is a well defined, unambiguous, and safe choice even for geometries far from the equilibrium. Unfortunately, this often leads to too large ansätze and therefore some simplification is necessary. In the case of ozone, one intends to use the CAS to generate an expansion space by applying single and double excitations to all determinants of the CAS for the treatment of static and dynamic electron correlation in a single wave function model. This can be done via a multireference CI Singles and Doubles (MR-CISD) ansatz in the MRCC program. In the case of a DBOC calculation, the implementation is limited to the use of single-reference HF-SCF orbitals, which represents, beyond the limited size of the CAS, one further compromise to be made. Nevertheless, it is reasonable to assume that the electron correlation effects in DBOC can still be well recovered, provided that the CAS is constructed properly.
For ozone, performing MR-CISD DBOC calculations with a fullvalence CAS (18,12) space is not possible, at least not in a reasonably sized basis set. In fact, even a CAS(6,6) space is computationally too demanding if, as discussed in the previous chapter, a triple-ζ quality basis set is to be employed.
Such challenges are well known in the field of high-accuracy ab initio thermochemistry [65][66][67][68][69][70][71][72][73][74] , where sub-kJ/mol accuracy is successfully achieved by considering the effect of higher and higher excitations in the electron correlation treatment via energy corrections evaluated in smaller basis sets. This approach, if carried out in the systematic hierarchy of the correlation consistent basis sets of Dunning and co-workers [75][76][77] , was found to work very well not just for contributions to the total energy, but also to structural and spectroscopic parameters 69,74,78 . Going this path with the DBOC surface for ozone, we define the function which thus accounts for the errors due to the incompleteness of the CAS(4,4) space with respect to the CAS(6,6) one via a correction term evaluated in the smaller cc-pVDZ basis set. The DBOC surface defined by eqn 2 is simply referred to as DBOC (6,6) in the subsequent chapters.

Computational Details
DBOC surfaces described by eqn 2 , as well as at the simple CAS(4,4)/cc-pVTZ levels were obtained in a set of points corresponding to the sparser grid of the high-accuracy PES from ref. 12. Altogether 329 points in the Born-Oppenheimer relative energy range of V BO < 7000 cm −1 , corresponding to structural parameters of 2.2 a. u. ≤ r 1 ≤ 2.6 a. u., 2.2 a. u. ≤ r 2 ≤ 3.0 a. u. and 100.0 • ≤ OOO ≤ 135 • were considered. The calculations were performed with the combination of the CFOUR 53 and MRCC 54 program codes at the MR-CISD level. The CAS(6,6) reference space consisted of the three highest-lying occupied and three lowest-lying unoccupied molecular orbitals of the 2p space (the space of orbitals formed by the 2p atomic orbitals of the oxygen atoms), while in the CAS(4,4) space only the highest-and lowestlying two occupied and unoccupied orbitals were included, respectively. Due to the lack of core polarization functions in the cc-pVXZ basis set series, all quantities had to be evaluated with the core electrons excluded from the correlation treatment. All calculations were carried out on the ROMEO 2018 supercomputer facility of the University of Reims Champagne-Ardenne. 79

Fitting the analytic models to ab initio DBOC values
The summary of calculated ab initio points representing DBOC values at a grid of nuclear geometries for the three symmetric isotopologues of ozone is given in Table 1. Here and below we use the standard convention for abbreviated notations of isotopic species (  internal bond length (r 1 , r 2 ) and apex angle (α) coordinates. A graphical overview of the DBOC corrections on the full grid is shown in Fig. 1 for the principal ozone isotopologue 16 O 3 . Differently colored groups of points correspond to two-dimensional cuts of the full PES for fixed values of the apex angle ranging from 100 • to 135 • . Several analytical ∆V DBOC (r 1 , r 2 , α) models involving polynomial, Gaussian, exponential, inverse hyperbolic functions and their combinations were tested to fit ab initio DBOC values on the full grid of nuclear displacements. At large nuclear displacements from the equilibrium, erratic deviations were obtained for all models. This erratic behavior of the DBOC PES with (ab initio -fit) discrepancies larger than several wavenumbers generally occurred at the energy range above the half of the dissociation energy threshold (D 0 / 2) , where the experimental value 80,81 of D 0 is about 8560 cm −1 ( as sited in 9). Also, for large nuclear displacements corresponding to V BO > D 0 / 2 we have experienced some convergence issues of ab initio DBOC calculations, probably to be attributed to the inappropriateness of single-reference HF-SCF orbitals in these structures.
The most accurate to date ab initio BO PES constructed by Tyuterev at al 12 provided vibrational levels calculations with the root mean square deviation of calculated and observed values (Observed-Calculated) of about 0.5 cm −1 for low energy vibrations and about 1 cm −1 up to 93 % of the dissociation threshold D 0 . The aim of this work being to improve this accuracy, we focus on the careful investigation of the DBOC corrections at various levels of the theory, at the BO PES range up to 4300 cm −1 that roughly corresponds to D 0 / 2. To this end, our grid of points were chosen significantly denser at the bottom of the potential well near the C 2v equilibrium (open geometrical configuration of the ozone molecule). Panels (B) and (C) of Table 1 give the corresponding information for the grids below 3000 cm −1 and below 4300 cm −1 , respectively. In the weighted fit we assigned the largest weights to the points of the (B) panel, with lower weights for the points presented in the (C) panel, whereas erratic points beyond D 0 / 2 were excluded from the final fit. Different analytical models gave us quite similar corrections to vibrational levels in this energy range. The final fits were performed with the following analytical form: where r e and α e are, respectively, the bond length and angle at the C 2v equilibrium structure, while a, b, c, C and d i jk are fitted parameters, and n is the order of the expansion. This functional form accounts for the C 2v symmetry of the PES for the considered isotopic species. The schematic behavior of the PES is described by the Gaussian function of r 1 , r 2 and cos α, and the finer tuning of the surface is done by means of the symmetrized polynomial part. The asymptotic behavior of the correction is controlled by the "global" Gaussian damping function which ensures that the correction approaches zero at sufficiently large values of r and α. The resulting isotope-dependent PESs then were constructed for each isotopologue separately by adding the fitted DBOC surface to the Born-Oppenheimer BO PES 12 in the form Here (i) denotes one of the 666, 777, and 888 isotopologues. The fitted models were coded in Python and can be found in the Supplementary Information. The residuals of the fit are collected in Fig. 2 for the principal isotopologue and in Table 2 for all three symmetric isotopologues. The deviations of the fit of the analytical ∆V DBOC PES (eqn 3) to ∆E DBOC values lie mostly within 0.2 cm −1 for the configuration energies below 2500 cm −1 , and within the 1 cm −1 corridor for the configuration with energies in 2500-4300 cm −1 range, as shown on the lower panel of Fig. 2.

Ab initio vibrational energy levels of 16 O 3 with DBOC corrections
Global variational methods permit nowadays converging vibrational basis set calculations to the precision of 0.01 or 0.001 cm −1 for small and medium sized molecules 12,24,[61][62][63][64][82][83][84] , at least in the energy range up to the half of the dissociation threshold. This level of precision is necessary for this study to clearly distinguish the contributions beyond the BOA because best published ab ini-tio BO ozone calculations 12,24,85 were already quite accurate with the average (Observed-Calculated) discrepancies of 0.5 cm −1 to 1 cm −1 . As often recommended, atomic masses were used for calculations of vibrational levels in order to partly account for nonadiabatic contributions. Here vibrational levels were computed in this range using the same variational technique as described in our previous work, 12 the only difference being the inclusion of the mass-dependent DBOC correction to the BO PES according to eqn 4. We consider the effect of vibrational DBOC corrections at two levels of the theory: using complete active spaces CAS (4,4) and CAS(6,6) as described in Section 2. The corresponding results for vibrational levels with the V BO + ∆V DBOC PESs will be denoted as DBOC (4,4) and DBOC(6,6) for the sake of brevity. The vibrational levels of the ozone molecule are usually assigned in terms of normal mode quantum numbers (v 1 ,v 2 ,v 3 ) where v 1 stands for symmetric vibration, v 2 for the bending mode and v 3 for the anti-symmetric vibration. Due to molecular symmetry and approximate coincidence of harmonic frequencies (ω 1 ≈ ω 3 ) they are organized in series of so called "stretching polyads" 3,22,25 . These polyads, denoted as P s,b , are defined as sets of nearby vibrational states where the ∪ symbol stands for a union over v 1 ,v 3 under the condition that v 1 + v 3 = s, which is the total quantum number of stretching vibrations, and b = v 2 is the bending quantum number. The intra-polyad couplings among vibrational modes in ozone are dominated by Darling-Dennison resonances and rotation-vibration couplings by Coriolis resonances 3,25 .
It is instructive to investigate the DBOC contributions to individual vibrational levels, as well as to the general polyad-bypolyad picture. An improvement of ab initio calculations for vibrational band centers with respect to observations is clearly seen in Tables 3 and 4 where the observed minus calculated discrepancies are collected together with the root-means-squares (RMS) and mean deviations for various types of polyads. Let us remind that all comparisons in this work are given up to energies about V BO =D 0 / 2 because of limitations in the construction of accurate DBOC surfaces described in Section 3. As it is shown in Table 4 and in Fig. 3, the inclusion of DBOC both at the CAS(4,4) and CAS(6,6) levels yield significant improvements with respect to the the most accurate available BO calculations 12 , with quite similar results for both active spaces.
The improvement is particularly pronounced for the bending states as follows from the first part of Table 3 and Fig. 5: the (Observed-Calculated) errors decrease by an order of magnitude when accounting for the DBOC(4,4) or DBOC(6,6) terms in the PES. For the (010), (020) and (030) bending vibrational states the BO+DBOC(6,6) surface provides nearly experimental accuracy as shown in Fig. 5. Note, that higher very weak bending bands 4ν 2 and 5ν 2 have not been directly observed. The corresponding vibrational states (040) and (050) were considered in experimental spectra analyses as "dark" ones 3 and have been evaluated indirectly 86 form resonance perturbations of certain ro-vibrational transitions. The experimental accuracy for these "dark" states roughly corresponds to the number of digits given in the first column of Table 3. Consequently, a part of the (Observed-Calculated) deviations for (040) and (050) have to be attributed to the experimental uncertainty. The same comments apply to other vibrational states that have been considered as "dark" ones in experimental spectra analyses. This is the case of (400) and (140) in Table 3, for which we give only one decimal digit for the virtually "observed" values. In such cases, the (Observed-Calculated) deviations can be impacted by experimental uncertainties. For a general view on the energy contributions involving stretching bands we plot in Fig. 6 and 7 the mean and RMS Table 3 Comparison of vibrational levels computed from ab initio potential energy surfaces at various levels of theory with experimental data for the 16   (Observed-Calculated) deviations for the G s band systems where s = v 1 + v 3 is a total stretching quantum number. Here G s band system gathers all P s,b polyads with various bending quantum numbers b = v 2 using the same cut-off as in Tables 3, 4 : where the ∪ symbol stands for a union over all b = v 2 values included in Tables 3 and 4. The mean vibrational energy of the G s groups of levels denoted as E mean (G s ) is shown along the horizontal scales of Figs. 6 and 7, in cm −1 . calculation errors averaged over bending quantum numbers as a function of excitations of the stretching ( valence ) vibrational modes. Again, it is seen that the inclusion of the BO breakdown corrections globally improve the set of stretching vibration levels for the considered energy cut-off. The BO+DBOC(4,4) and BO+DBOC (6,6) show again very similar performance with the average (Observed-Calculated) discrepancies below 0.5 cm −1 .

Ab initio vibrational energy levels of 18 O 3 and 17 O 3 isotopologues with DBOC corrections
As the next step of this work, let us now examine the mass dependence of the DBOC corrections in vibrational levels for the homogeneous isotopic substitutions 18 Tables 5 and 6, as well as in Figs. 8-11. As for 18 O 3 , we obtain the most pronounced improvement for the bending states -by more than an order of magnitude -even though the observed series is much shorter, limited here to the ν 2 and 2ν 2 bands. The inclusion of DBOC corrections also permitted significantly more accurate results for mixed bendstretch polyads: the RMS (Observed-Calculated) deviation drops down from 0.93 cm −1 to 0.39 cm −1 , while the mean (Observed-Calculated) deviation is reduced from 0.77 cm −1 to no more than 0.04 cm −1 . For pure stretching polyads the improvement is less pronounced, possibly because the errors obtained with the BO PES 12 in the ν 1 band and in the combination series ν 1 + nν 3 series were already very small (between 0.02 cm −1 to 0.2 cm −1 ). This applies for 16 O 3 , 18 O 3 and 17 O 3 species as well. However, the general trend in improvement of the total (Observed-Calculated) statistics is clearly confirmed by the reported comparisons for all three isotopologues.  Table 7.

Discussion and conclusions
In this work, we computed mass-dependent diagonal Born-Oppenheimer corrections to the ab initio electronic ground state potential energy surface for the main 16 O 3 isotopologue of the ozone molecule, as well as for homogeneous isotopic substitutions 17 O 3 and 18 O 3 , which had not been investigated previously.
The goal was to improve the accuracy of the PES with respect to currently best results for vibrational states obtained in the BO approximation 12 which was in average about 0.5 cm −1 for low vibrations and about 1 cm −1 for the entire set of observed bands. The other objective was to investigate the mass-dependence of the PESs in the 16 O 3 , 17 O 3 , and 18 O 3 isotopic series and the corresponding impact on vibrational states. As the ozone molecule is known to be a strongly multireference system, we employed to this end an MR-CISD ansatz in the MRCC program with CAS (4,4) and CAS (6,6) actives spaces. The DBOC contributions to electronic energies were then fitted to an analytical form in terms of nuclear displacements from the equilibrium. Finally, the vibrational band centers were computed using nuclear motion variational method using the BO+DBOC PESs for the three species     the accuracy with respect to the best BO calculations. For the low vibration range the BO calculations of ref. 12 using atomic masses were generally slightly underestimated with positive (Observed-Calculated) deviations. The DBOC(4,4) and DBOC(6,6) corrections modify this in the right direction, providing a remarkable reduction of the calculation errors as is clearly seen in Tables 3-4 and on Figs. 5-7. The account of DBOC produces the most pronounced improvement for bending states and for mixed bendstretch polyads. In the case of (010), (020) and (030) states the errors are reduced by more than order of magnitude: the (Observed-Calculated) deviation for BO+DBOC(6,6) are below 0.02 cm −1 (Fig. 5 and Table 3), approaching nearly experimental accuracy. Similar improvement occurs for the bending states of 18 O 3 (Table 5). For combination states with three vibrational quanta, the DBOC corrections bring ab initio calculations to the mean errors below 0.2 cm −1 as shown in Tables 3-6 and Figs. 3-11. To our knowledge, this level of ab initio accuracy has not been obtained so far for molecules with such a complex multireference electronic structure as ozone. Note that DBOC (4,4) appears to be competitive with DBOC(6,6) and provides even slightly better results for some polyad series. Possibly, the convergence with respect to the size of the CAS in the considered energy range was approached. Let us remind, however, that full variational electronic energy calculation with the triple-ζ quality basis set were only performed with the the CAS(4,4) space. CAS(6,6) calculations were only possible using an approximate incremental scheme (eqn 2) using a double-ζ basis set, a novel approach clearly justified by the results. It is curious to note in Figs. 3 and 6 that the mean (Observed-Calculated) deviations would tend to almost zero values if one takes the average of DBOC(4,4) and DBOC(6,6) -green and red lines in Figs. 3,6 -which, although possibly no more than a fortunate coincidence, can be interpreted as a hint that the incremental formulation defined in eqn 2 might be overestimating the effect of the larger CAS to a certain extent. Anyway, in terms of the mean (Observed-Calculated) deviations, the improvement is quite striking. In total, they drop down from 0.69 cm −1 (BO) to 0.05 cm −1 (DBOC(4,4)) or to -0.12 cm −1 (DBOC(6,6)) for the considered range including the set of 16 O 3 vibrations in Table 3 with five bending and four stretching quanta. For 18 O 3 , the mean (Observed-Calculated) deviations drop down from 0.55 cm −1 (BO) to -0.09 cm −1 (DBOC(6,6)). This means that DBOC corrections almost remove the systematic BO errors in vibrational levels below E BO = D 0 / 2. The remaining (Observed-Calculated) scattering naturally increases with energy that could be attributed both to BO and to DBOC uncertainties in the analytical models, as well as to individual vibrational modes. However, the DBOC corrections generally reduce the total RMS (Observed-Calculated) deviations by a factor of two or so: from 0.84 cm −1 to 0.41 cm −1 for 16 O 3 , from 0.73 cm −1 to 0.41 cm −1 for 18 O 3 and from 0.49 cm −1 to 0.26 cm −1 for 17 O 3 . The present results, including full sets of calculated vibrational states provided in the Supplementary Information with (Observed-Calculated) statistics, can be used for improving spectroscopic effective models for analyses of observations. Despite a significant progress 3,26 , some ozone bands below E BO = D 0 / 2 have not yet been fully analyzed. For complicated coupled band systems belonging to overlapping polyads 92 , ab initio values are often used for constraining parameters which are missing from direct measurements. Typically, "dark" states are fixed to theoretically predicted values. Their accuracy has thus a considerable impact on spectral perturbations induced by accidental rovibrational Coriolis or Darling-Dennison resonances 3,25,26,85 . This is particularly important for rare species like 18 O 3 and 17 O 3 , were much fewer spectral assignments are yet available . As explained in sections 2 and 3, reliable DBOC calculations with the targeted accuracy were possible to conduct up to about the half of the dissociation threshold -at least using the currently feasible ab initio ansatz with accessible computational facilities. Above this energy corresponding to ≈ 4300 cm −1 , ab initio calculations were not well converged for certain nuclear geometries and / or the fits to analytic models yielded erratic series with outliers much larger than 1 cm −1 despite the use of several trial functions. Although this could be due to an insufficiently large active space or a limited cardinal number for the atomic basis set, we consider the limitation to the use of single reference Hartree-Fock Self-Consistent-Field orbitals in the DBOC calculations as the major obstacle to achieving higher accuracy. Note that the most accurate to date BO PES 12 , which we used in this work, had been constructed using CAS (18,12) with Complete Active Space Self-Consistent Field (CASSCF) orbitals in the X=5,6 → ∞ extrapolated complete basis set limit. In addition, it accounted for a supplementary correction 12 due to a simultaneous optimization of 13 excited electronic states considered by Dawes et al 10 , which had a significant impact on shape of the PES 10,11,14 and the dynamics 35,37,38 in the range near the transition state towards the dissociation. Another challenge is to account for non-adiabatic interactions among 27 electronic states tending to the same dissociation limit 39,93 and for possible effects of the topological phase 8 related to conical intersections 40 . However, the reported results at the present state of art already provide an encouraging accuracy validation for the multireference methods of the ab initio theory, at least in the range up to four vibrational quanta. In future works, we plan to address the above issues, as well as to extend DBOC calculations to nonhomogeneous isotopic species of ozone.

Conflicts of interest
The authors declare no conflicts of interest.
download file view on ChemRxiv paper.pdf (2.07 MiB)