Testing the limitations of harmonic approximation in the determination of Raman intensities

Raman intensities in molecular spectra are usually computed within double harmonic approximation. This procedure relies on treating a vibrating molecule as a collection of harmonic oscillators and on the assumption that polarisability tensor invariants display linear variations around the molecular equilibrium geometry. This methodology, originally formulated by Placzek, constitutes the theoretical foundation for computing Raman intensities in standard quantum chemistry packages. However, the two assumptions underlying double harmonic approximation have not been sufficiently tested. In this work, we employed exact anharmonic ro-vibrational wave functions and distance-dependent polarisability invariants together with their harmonic approximants to investigate the discrepancies in Raman intensities of the fundamental transitions in 12 diatomic molecules, caused by double harmonic approximation. We found that: (i) the errors in total Raman intensities were between −8.2% and +9.5%, (ii) the largest discrepancy was observed for F2, where the polarisability invariants could not be adequately modelled by their linear approximants, and (iii) quantum chemical methods fail to predict reliable polarisability invariants at non-equilibrium molecular geometries; the associated errors in Raman intensities are huge and completely overshadow the shortcomings of double harmonic approximation. We communicate here an urgent need for developing accurate methods capable of computing reliable polarisabilities also at distorted geometries. GRAPHICAL ABSTRACT


Introduction
Raman spectroscopy is a well established branch of vibrational spectroscopy. Theoretical determination of Raman transition frequencies and the corresponding intensities comprises an important part of spectroscopic analysis. Theoretical foundations for this procedure were laid by Placzek [1] who used the polarisability tensor and vibrational wave functions as building blocks of his theory. In his work, the first-order induced dipole P ρ along a Cartesian axes ρ, is expressed as a function of the polarisability Q, α ρσ ≡ α ρσ (Q). Polarizability is then expanded at the equilibrium molecular geometry as a Taylor series in Q to obtain α ρσ = (α ρσ ) 0 + k δα ρσ δQ k 0 Q k + 1 2 l k δ 2 α ρσ δQ k δQ l 0 Q k Q l + · · · (3) Substituting the above expansion for polarisability truncated at the first-order, and using harmonic oscillator wave functions for ψ in Equation (2) while working under the assumption of electrical harmonicity gives (4) which consists of two terms. 1 Using selection rules derived from harmonic-oscillator functions (see Section S1 in the supplementary material for more details), the zeroth-order term corresponds to Rayleigh scattering, while the latter term, containing the first-derivatives of polarisability corresponds to Raman scattering.
The above discussion based on the harmonic approximation constitutes the primary approach to the computation of Raman intensities [2]. 2 Both formalisms for normal coordinate analysis and the resulting intensity analysis rest on the harmonic approximation. Consequently, almost all quantum chemistry programs available at present employ this double harmonic approximation for the computation of vibrational intensities, with truncation of polarisability invariants at the first order. While a great deal of work has been done regarding the discrepancy of harmonic transition frequencies with experimental results [3][4][5][6][7][8][9][10][11], a similar discussion pertinent to Raman intensities is limited [12,13]. The question to be answered is: What is the error associated with the assumption of the electrical harmonicity? In other words, how accurate are the Raman intensities computed under the harmonic approximation?
There are a few challenges to overcome in order to reliably answer these questions. The first one is to obtain accurate state specific wave functions for the studied molecules, and the second is to obtain reliable polarisabilities expressed as a function of the vibrational coordinates. These issues limit the choice of the possible molecules which can be studied in this context. For diatomic molecules, the exact state-specific vibrational wave functions can be readily obtained beyond the harmonic approximation. In contrast, for larger molecules this is significantly harder. In computation of ab-initio polarisabilities, the accuracy of the results depends to a large extent on the choice of the quantumchemical technique used for the calculations and on the choice of the basis sets used to construct the electronic wave function. Thus, it becomes computationally challenging to reliably perform tests on the harmonic approximation on larger molecular systems, where the quantum chemical calculations become expensive or even infeasible. Thus, smaller molecules are preferred for such a study. In the case of diatomic molecules, the polarisability can be naturally expressed as a function of the internuclear distance. In view of these considerations, we limit ourselves in the current study to a series of diatomic molecules containing from two to 18 electrons.
For diatomic molecules, several experimental and theoretical studies have been reported in the literature which indirectly expose the errors originating from the double harmonic approximation. Herman and Wallis discussed the centrifugal distortion arising from the vibration-rotation coupling during accurate determination of infrared intensities [14]. Similar study for the Raman intensities was performed by James and Klemperer [15]. In these works, the diatomic molecule under study was assumed to behave as a harmonic oscillator, and the centrifugal term J(J + 1)/2μr 2 was approximated as a 2-3 term expansion over the reduced internuclear distance, ξ , defined as ξ = (r − r e )/r e . Polarizability anisotropy γ = α || − α ⊥ , an invariant of the polarisability tensor α, governing the intensities of the O-and S-branches was approximated as a Taylor series expansion over ξ γ (ξ) = γ 0 + γ 1 ξ + 1 2 γ 2 ξ 2 + · · · where γ 1 = ∂γ ∂r and γ 2 = ∂ 2 γ ∂r 2 , while retaining up to three terms in the expansion. Using perturbation theory or numerical solutions, the correction factors (usually referred to as the Herman-Wallis factors) were approximated as functions of derivatives of polarisability anisotropy (γ n ), rotational constant (B e ), vibrational constant (ω e ) and the rotational quantum number (J). Alternatively, using known values of molecular constants and the determined Herman-Wallis factors, the values of derivatives of polarisability anisotropy for equilibrium internuclear distance were obtained by analysing intensities from acquired Raman spectra. The first derivative of polarisability anisotropy (γ 1 ) for H 2 and D 2 were determined by Asawaroengchai and Rosenblatt [16]. Hamaguchi and coworkers employed a simplified expansion [17,18] given by and reported the numerical value of (γ 1 /γ 0 ) r e and (γ 2 /γ 0 ) r e for the same molecules. In this approach, the calculated Raman intensities were fitted using non-linear least squares regression to experimental intensities for O0-and S0-branch lines of H 2 and D 2 , obtaining the values of (γ 1 /γ 0 ) r e and (γ 2 /γ 0 ) r e as fit parameters along with the intensity response curve of the spectrometer. Further, a similar analysis performed for N 2 revealed that the inclusion of the second derivatives in Equations (5) and (6) is less important than for H 2 [17][18][19].
In our earlier work, we have investigated the effects of centrifugal distortion to the Raman intensities of purely rotational transitions in H 2 , HD and D 2 , and examined the importance of including higher order derivatives of polarisability invariants in the process of determination of the Raman intensities [20]. This work revealed that Taylor series expansion up to the second order yields an excellent representation of both the invariants around r e . We start our analysis by performing accurate computations of Raman intensities determined from accurate, anharmonic ro-vibrational wave functions of diatomic molecules and the corresponding geometry-dependent molecular polarisability invariants computed over the varying internuclear distance using standard ab initio methods coupled with response theory [21][22][23][24][25][26][27]. Accurate ro-vibrational wave functions of diatomic molecules are obtained by numerical solution to the radial nuclear problem given by Equation (7). These quantities allow us to evaluate numerically the integral in Equation (2) in order to compute the transition matrix element for a specific transition. In the next step, we construct approximations to the distance-dependent polarisability invariants by constructing their Taylor expansions at the internuclear distance corresponding to the ground rovibrational state in Equation (3) and truncating them at some low order. Such simplified polarisability invariants representations are subsequently used in the spirit of original Placzek approach to determine the Raman intensities using the integral in Equation (2). Further, we analyse the role of an accurate representation of the rovibrational wave functions for the computation of the integral in Equation (2), comparing the Raman intensities computed with accurate anharmonic ro-vibrational wave functions with those computed using the corresponding harmonic ro-vibrational wave functions.
The current study relies to a large degree on the knowledge of accurate potential energy curves used to determine numerically the vibrational wave functions of diatomic molecules. In particular, the underlying potential energy curve for the molecular hydrogen, being one of the pillars of the current analysis, is known to high accuracy owing to the extensive work of the late Professor Lutosław Wolniewicz. He devoted many decades of his life for producing more and more accurate and physically exact potential energy curves for H 2 and its isotopologues. With this publication we would like to pay a tribute to the memory of Professor Wolniewicz, without whose magnificent work the present analysis could not be adequately performed.

Overview
Two components are required for accurate evaluation of the integral in Equation (2) for diatomic molecules: (a) polarisability and (b) ro-vibrational wave functions, both expressed as a function of the internuclear distance. In this section, we define the necessary terms relevant to the present analysis.
Polarizability tensor defined in the Cartesian coordinates consists of two non-vanishing principal components, α ⊥ , defined as the component perpendicular to the molecular axis, and α , defined as the component parallel to the molecular axis. For freely rotating molecules, two rotational invariants of polarisability are employed instead: mean polarisabilityᾱ = (2α ⊥ + α )/3 and polarisability anisotropy γ = α − α ⊥ .
Due to the simple structure of diatomic molecules, their accurate ro-vibrational wave functions can be expressed as functions of a single variable, the internuclear distance (r). The ro-vibrational wave functions are determined by numerical solution to the radial nuclear problem given by Equation (7). Each ro-vibrational function obtained in this way is indexed with the vibrational quantum number v and a rotational quantum number J, and-for our convenience pertaining to the numerical integration procedure used here-also by the span of the internuclear distances over which the function is non-negligible.
In order to generate reliable reference data for testing the double harmonic approximation usually used for computing Raman intensities, 3 we start by evaluating the integral in Equation (2) in the most accurate form, i.e. with anharmonic ro-vibrational wave functions ψ anhrm and exact polarisability invariants (ᾱ exact or γ exact ). Both quantities are computed over a fine grid of internuclear distances and represented as cubic splines. Raman intensity for the fundamental transition is then computed using the relevant transition matrix elements. This quantity constitutes our reference against which the approximate forms of the same quantity are compared.
In the next step, we construct approximate forms of the quantities needed to evaluate the integral in Equation (2), namely, sets of harmonic wave functions ψ harmn determined from harmonic potentials and Taylor series expansions of polarisability invariants at the equilibrium internuclear distance, denoted asᾱ (n) or γ (n) , where where n refers to the Taylor series expansions truncated at the n th order. The approximate quantities are then used to compute the relevant matrix elements, and finally the corresponding Raman intensities.
We conclude our analysis by performing a detailed comparison between the exact and approximate sets of Raman intensities of the studied diatomic molecules and present a discussion of accuracy of the studied approximations.

Polarizability calculations
Static (frequency-independent) components of the polarisability tensor were computed over a fine grid of internuclear distances for all the studied molecules using various ab initio techniques. Specific details for each molecule are given below.
Distance-dependent static polarisabilities of H 2 , HD and D 2 were adapted from our earlier work [28]. In brief, the polarisability was computed using CCSD technique [29,30] (equivalent in this case to FCI calculations) using a composite basis set comprising of hydrogen's atomic aug-mcc-pV6Z basis set of Mielke et al. [31] downloaded from the EMSL basis set database [32,33] and custom designed bond-functions. These calculations were performed using DALTON quantum chemistry package [30]. In the case of H 2 , HD and D 2 , a large discrepancy in the values of mean polarisability and its anisotropy were observed when using standard atomic-basis in comparison with the results computed using explicitly correlated wave functions by Rychlewski [21][22][23]. The application of standard AO basis in conjunction with additional bondfunctions placed at five equidistant positions along the internuclear axis resulted in significantly improved polarisabilities, comparable in quality to the explicitly correlated results of Rychlewski. The additional bond functions comprised of even-tempered Gaussian primitives (8s6p), constructed using the expression αβ 1−k with α = 1.8518519 and β = 3 for the s-functions and with α = 2.520 and β = 3 for the p-functions. The detailed development of this composite basis set was described in our earlier work [28]. The previously presented results filled up gaps in the available data relevant to the spectroscopy of molecular hydrogen, and resolved several irregularities in the polarisability datasets reported by Rychlewski [22].
The polarisabilities for the remaining molecules studied here (HF, HCl, CO, N 2 and F 2 ) were initially computed using response theory with the coupled cluster with single, double and perturbative triple excitations (CCSD(T)) wave functions [34,35] and the aug-cc-pVQZ basis sets [36,37]. The cfour quantum chemistry package was used for these computations [38,39]. Keeping in mind the development for the molecular hydrogen described above, it was natural to expect that additional bond functions placed along the internuclear axis would systematically improve the computed polarisabilities also for HF, HCl, CO, N 2 and F 2 . Unfortunately, this was not the case. The gradually developed bond-functions (for more details, see Section S2 of the supplementary material) allowed us to systematically improve the total energy values for these molecules, but at the same time the observed changes in the polarisability invariants displayed rather unsystematic and at times unexpected convergence patterns with no apparent trends. Consequently, we gave up the idea of using composite basis sets and decided to use the CCSD(T)/aug-cc-pVQZ polarisabilities in our analysis.

Wavefunctions
The ro-vibrational wave functions ψ v,J ≡ ψ v,J (r) for the studied diatomic molecules were obtained by numerically solving the radial nuclear equation where the potential energy curve V(r) was adapted from earlier works or determined in a way described below using reported sets of experimental transition frequencies.
The potential energy curves for molecular hydrogen and its isotopologues were adapted from the work of Wolniewicz [57] who reported the Born-Oppenheimer potential together with the scaled adiabatic, relativistic, and radiative corrections to it. Each of these components was interpolated to a finer grid spanning the interval from 0.2 to 12.0 a.u. and combined together in order to obtain the final potential energy curve. These curves were given an analytic, third-order spline representations, which was subsequently used in numerical solution to the radial nuclear Schrödinger problem. For details, see our previous work [28].
The usual way of determining accurate potentials V(r) for solving radial nuclear Schrödinger problem is based on quantum chemical calculations. Here, however, we resorted to a different solution, which can be referred to as an inverse Schrödinger problem. To this end, we employed a set of accurate experimental ro-vibrational levels and transition frequencies to determine the best potential energy curve that reproduces these levels and frequencies via the solution to Equation (7). There are many possible choices of the analytical form for the potential energy function that can be used in this context. Here, for the HF, DF, HCl, DCl, N 2 , 14 N 15 N, and F 2 molecules, we chose the modified Morse function [58,59] with five fitting parameters α, c 0 , c 1 , c 2 , and r e and for CO and 13 C 16 O, we used the Rydberg function [60-62] with b = k e D e and with three fitting parameters k e , D e , and r e . The optimal values of the parameters for each potential energy function were obtained by performing a non-linear least squares optimisation with experimental ro-vibrational transition frequencies and energy levels covering v = 0-3 and J = 0-4 used as reference. The optimised parameters are tabulated in Table 1. A detailed description of the optimisation procedure, along with the list of experimental ro-vibrational levels used in the optimisation and the corresponding values obtained by solving the ro-vibrational Schrödinger equation with the fitted potentials, are given in Section S3 of supplementary material.
A similar approach was used to determine the corresponding harmonic wave functions ψ harmn . The only difference to the just discussed scheme used for finding accurate ro-vibrational wave functions ψ anharm concerned using a set of harmonic energy levels for each molecule as a reference and an analytical potential expression V(r) = k(r − r e ) 2 with r e corresponding to the equilibrium internuclear distance with the force constant k optimised to reproduce the harmonic vibrational frequency ω e of each system. The values of r e and ω e were taken from the compilation of Herzberg and Huber [63]. The values of these constants together with the reproduced harmonic levels are given in Section S4 of supplementary material.
The solution to the radial nuclear Schrödinger equation in Equation (7) was constructed numerically using the collocation method [64,65] with the five-point stencil representations of the derivative operators. The Hamiltonian matrix was diagonalised for each rotational state J separately. A detailed process of constructing the Hamiltonian matrix and its diagonalisation is described in detail in Section S4 of the supplementary material. The resulting ro-vibrational wave functions ψ v,J were determined only on the collocation grid and their analytical representation was subsequently obtained as a third-order spline. The relevant scheme for solution of Equation (7) was implemented in Python using NumPy [66,67].
The spline representations [68] of the computed polarisability invariants =ᾱ, γ and the resulting ro-vibrational wave functions ψ v,J were subsequently used in an adaptive Gauss-Kronrod-Patterson quadrature [69][70][71] to compute numerically the integral in Equation (10) yielding the required transition matrix elements The lower and upper integration limits, r min and r max , respectively, were selected in such a way that both the involved ro-vibrational wave functions ψ v,J and ψ v ,J were negligible (|ψ v,J |, |ψ v ,J | < 10 −6 ) beyond the interval [r min , r max ]. The integral in Equation (10) can be used for computing the normalisation of each ro-vibrational wave function if we set = 1. We use this condition to determine Table 1. Details regarding the potential energy curve for the studied molecules in this work. For hydrogen and its isotopologues, we adapt the potential energy curve from the work of Wolniewicz (Ref. [57]) after cubic-spline interpolation (Ref. [28]). For the other molecules, the optimised parameters of the analytical potential functions used are listed (see text for more details). Ref. [28,57] r e = 2.13081292 The parameter c 0 was sufficient to reproduce the observed energy levels in this case, and hence, c 1 and c 2 were set to zero.
the normalisation constants for each of the used here ro-vibrational wave functions, so in the following equations and discussions all the normalisation constants will be automatically equal to one and omitted. For interested reader, we mention that all those wave functions are available from GitHub [72] repository; for details, see Section 6.

Truncated expansions of polarisability invariants
Truncated Taylor series approximations of polarisability invariants at the internuclear distance r 0 = ψ v=0,J=0 | r |ψ v=0,J=0 were constructed as follows: • The values of =ᾱ, γ computed on the original grid were fitted over the interval [r min , r max ] using a nine-degree polynomial. We used typically 50-60 grid points for this step. The resulting polynomial was represented in a Taylor-series-like form where (n) r 0 can be interpreted as the n th order derivative of at r 0 . The fit reproduced exact faithfully with typical deviations of around 1×10 −4 ( ≤ 0.1%).
• Only exact and (1) were used to compute the matrix elements of the polarisability invariants in conjunction with the exact anharmonic wave functions ψ anharm or their harmonic approximations ψ harmn . Squares of these matrix elements are tabulated in Tables 2, 3, and 4.

Results
The present analysis relies to a large degree on accurate, distance-dependent polarisabilities and ro-vibrational wave functions as the necessary ingredients. Therefore, we start our discussion by presenting the computed perpendicular and parallel components of the polarisability tensor for each molecule together with the corresponding rotational invariants. These results are followed by a short discussion of the ro-vibrational wave functions and a compilation of the computed matrix elements and the resulting Raman intensities.

Polarizabilities
The distance-dependent polarisability components (α ⊥ and α ) and the associated invariants (ᾱ and γ ) for molecular hydrogen are shown in Figure 1. These results, adapted from our previous work [28], were found to be in good agreement with the analogous results computed Table 2. Comparison of the matrix elements of polarisability invariants and the total Raman intensity a for the fundamental vibrational transition |v = 1, J = 0 ← |v = 0, J = 0 for molecular hydrogen and its isotopologues computed using CCSD methodology, and with exact anharmonic wave functions (ψ anhrm ) or their harmonic counterparts (ψ harmn ), along with exact polarisability invariants (ᾱ exact and γ exact ) or their linear approximations (ᾱ (1) and γ (1) ).
4%) a Corresponding to linearly polarised excitation and parallelly and perpendicularly polarised detection scheme, usually denoted in literature as (incident , detection +⊥ ). Table 3. Comparison of matrix elements of mean polarisability for the fundamental vibrational transition |v = 1, J = 0 ← |v = 0, J = 0 for nine diatomic molecules, computed at various levels of theory (HF, CASSCF, DFT and CCSD(T)) and approximations: with exact anharmonic ro-vibrational wave functions (ψ anhrm ) or their harmonic counterparts (ψ harmn ), and with exact mean polarisability (ᾱ exact ) or its linear approximation (ᾱ (1) ). using explicitly correlated wave functions by Rychlewski [22]; for more details, see [28]. The computed polarisability tensor components of molecular hydrogen at large internuclear distances converge properly to the correct atomic limit equal to twice of the atomic polarisability of hydrogen [73].
For the other studied here molecules (HF, HCl, CO, N 2 , and F 2 ) the analogous results are shown in Figures 2   and 3. Since the electronic wave functions were computed under the Born-Oppenheimer approximation, the resulting polarisabilities are the same for pairs of isotopologues: HF and DF, HCl and DCl, CO and 13 C 16 O, and N 2 and 14 N 15 N. No accurate reference data is available for these molecules in literature, so to be on the safe side we used various quantum chemical methods to verify the correctness of the presented here results. In the Table 4. Comparison of matrix elements of polarisability anisotropy for the fundamental vibrational transition |v = 1, J = 0 ← |v = 0, J = 0 for nine diatomic molecules, computed at various levels of theory (HF, CASSCF, DFT and CCSD(T)) and approximations: with exact anharmonic ro-vibrational wave functions (ψ anhrm ) or their harmonic counterparts (ψ harmn ), and with exact polarisability anisotropy (γ exact ) or its linear approximation (γ (1)   case of α ⊥ , the results obtained using HF, CASSCF, DFT, CCSD(T) and (at times) CCSDT show similar trends in the interval [r min , r max ], with the values increasing gradually with increasing internuclear distances. We conclude that the computed values of α ⊥ can be considered as reliable, but their accuracy cannot be expected to be larger than a few percent as substantial numerical discrepancies between curves obtained with various methods suggest.
In the case of α , the situation is much more serious. In principle, the curves computed with various ab initio methods for CO, HF, N 2 , and HCl show similar trends and resemble each other, except maybe for the internuclear distances close to r max , where the discrepancies between the curves are larger. However, for F 2 F 2, considerable deviations between the computed curves are observed inside almost the whole interval [r min , r max ], with particularly large discrepancies at larger internuclear distances. This should not be surprising as a similar observation was reported previously by Maroulis [74], who studied the distance-dependent polarisability and hyper-polarisabilities of F 2 F 2 using different ab initio techniques. From our perspective, however, of the planned analysis of Raman intensities, this is a rather unexpected and unforeseen obstacle, as it is difficult to say without further analysis which of these curves are correct. To resolve this issue we decided to look deeper into this problem by studying the distance dependence of α and α ⊥ at the full range of internuclear separations. The results of this investigation are presented in the next section.

Sanity check for the computed distance-dependent polarisabilities
The values of α and α ⊥ computed over larger range of internuclear distances (from 1.2 to 7 bohr) are shown in Figure 4 for all the studied here molecules except molecular hydrogen. At large distances, the values of both components of molecular polarisability should converge to the corresponding atomic limit equal to the sum of atomic polarisabilities. These limits, computed from accurate atomic data [73], are depicted in Figure 4 using purple triangles located at the right edge of each panel. Obviously, most of the used here methods are singlereference in nature and cannot describe adequately the process of molecular dissociation. Therefore, we do not expect that the curves computed with these methods (HF, DFT, and CCSD(T)) would converge to the correct limits. The question we try to address here is rather when these methods loose their applicability and start diverging.
Rather surprisingly, we discovered that α ⊥ computed with all the tested here ab initio methods converge-more or less accurately-to the expected atomic limit. In case of CCSD(T), the calculations for CO and N 2 cannot converge for internuclear separations larger than 4.6 bohr, but even so, the values of α ⊥ at 4.6 bohr are close to the corresponding atomic limits. The situation is completely different for α . Practically all the tested here singlereference ab initio methods diverge away from the atomic limits, but the trend of divergence is quite different for each method.
The only method which passes our sanity check is CASSCF. Again, this is not surprising, as CASSCF with valence active space is capable of producing correct wave functions over the whole range of internuclear separations, and consequently reproducing physically meaningful values of both components of the polarisability tensor, which converge to the correct atomic limit. Unfortunately, despite of quite reasonable distance dependence of α and α ⊥ computed with CASSCF, we proceed rather cautiously with accepting the CASSCF values in our analysis of Raman intensities, because the CASSCF polarisabilities seem to differ from more accurate CCSD(T) and DFT values at distances where all the methods are applicable. The best example can be inferred from Figure 2, where α ⊥ of F 2 computed with CASSCF and HF differ from the values obtained with other methods by 5-15%. We believe that reliable distance-dependent polarisabilities of many-electron systems-and consequently also reliable Raman intensities of such systems-can can computed only using multireference correlated methods, but the appropriate response theory tools are not yet available. We are forced to use in our further analysis the CASSCF values of the polarisability tensor invariants as the reference data even if we know that these values are not too accurate. A comparison of results obtained with other methods is presented later, in Section 3.4, where we use the CASSCF results as a reference and try to estimate the deviations in the total Raman intensities obtained with other ab initio techniques.
For some molecules, we observe spiky local features (α of N 2 at 3.5 bohr with the HF method and α ⊥ of F 2 at 6.0 bohr with the CASSCF method) or discontinuities (α and α ⊥ of CO at 4.1 bohr with the CCSD(T) method) in the computed polarisabilities in Figures 2-4. These originate from the shortcomings of response theory at stretched geometries. Fortunately, all of these spurious features are quite local and do not affect the polarisabilities in the region important for our analysis.

Ro-vibrational wave functions
Ro-vibrational wave functions ψ anharm for H 2 , HD and D 2 were obtained by a numerical solution of the radial nuclear equation with accurate potential energy functions determined by Wolniewicz [57]. The computed energy levels of ro-vibrational states are compared with accurate experimental and theoretical results in Section S4 of supplementary material. This comparison shows deviations within one wavenumber.
For the other molecules studied in this work, no potential energy curves have been reported in literature that would allow to achieve similar accuracy in the determination of ro-vibrational wave functions ψ anharm . Therefore, to circumvent this difficulty we decided to use an alternative solution that can be referred to as the inverse Schrödinger problem, in which accurate experimental ro-vibrational transition frequencies reported in literature were used to find an optimal shape of the potential energy surface allowing to reproduce these values with best fidelity; more details are given in Section 2.2.2. All the ro-vibrational energy levels computed with this procedure together with all the relevant transition frequencies were found to be within one wavenumber from the reference values. Detailed list of these energy levels and the corresponding transition frequencies is given in Section S4 of supplementary material together with their comparisons to experimental reference data.
The following conclusions can be drawn: (1) All the approximations, either pertaining to the simplifications in the ro-vibrational wave functions or to the simplifications in the polarisability tensor invariants, invariably lead to reducing of the values of the computed matrix elements. This eventually results in reduced total Raman intensities with respect to the exact values. This effect is more pronounced when the exact ro-vibrational wave function is replaced by its harmonic approximation (error of 6-8%). Approximating the exact polarisability by its linear form gives smaller error of about 1%. (2) The deviations from the exact results in the matrix elements of polarisability anisotropy are always a few times larger than analogous deviations in the matrix elements of mean polarisability, for all types of approximations. (3) The largest deviation (16%) in the computed matrix elements is observed for the polarisability anisotropy of H 2 computed by combining harmonic wave function ψ harmn with exact form of the invariant. Surprisingly, using linear approximation to polarisability anisotropy reduces this deviation by about half, owing to error cancellations. (4) The deviations in Raman intensities are comparable to those observed for the matrix elements of mean polarisability, due to the fact that the contribution from the polarisability anisotropy matrix elements-laden with substantially larger deviations-is diluted by a small prefactor of 7/45. (5) The double harmonic approximation, corresponding to the choice of ψ harmn andᾱ (1) and γ (1) in the integral in Equation (10), gives an overall error of around 2-3% in the computed Raman intensities. The deviation is slightly smaller for heavier isotopologues of molecular hydrogen.
We have analysed the data for molecular hydrogen separately above, because these results are not affected by possible large systematic errors in the determination of the polarisability invariants used to compute the Raman intensities. We believe that the results for molecular hydrogen are very close to the exact quantum mechanical limit. The results for other molecules discussed below use quite approximate values of the polarisability tensor components as discussed earlier in Section 2.2.1. The effect of the double harmonic approximation can be still approximately assessed using these CASSCF invariants, but the absolute values of the resulting Raman intensities should Table 5. Comparison of the total Raman intensity for the fundamental vibrational transition |v = 1, J = 0 ← |v = 0, J = 0 for nine diatomic molecules computed at various levels of theory (HF, CASSCF, DFT and CCSD(T)) and approximations: with exact anharmonic wave functions (ψ anhrm ) or their harmonic counterparts (ψ harmn ), along with exact polarisability invariants (ᾱ exact and γ exact ) or their linear approximations (ᾱ (1) and γ (1) (1) , γ = γ (1) 0 be treated with caution. We hope that more accurate means of determination of distance-dependent polarisability invariants will be available soon and these numerical results can be rectified in this way.
The results for HF, DF, HCl, DCl, CO, 13 C 16 O, N 2 , 14 N 15 N, and F 2 are presented in Table 3 (mean polarisability matrix elements), Table 4 (polarizability anisotropy matrix elements), and Table 5 (total Raman intensities). This time, the results are computed with three different approaches: (a) the exact approach with ψ anharm andᾱ exact and γ exact , (b) single harmonic approximation with ψ anharm andᾱ (1) and γ (1) , and (c) doubly harmonic approximation with ψ harmn and α (1) and γ (1) . The following discussion is based on the CASSCF results, but for completeness, Tables 3-5 contain data computed also with other methods (HF, DFT, and CCSD(T)). These results are discussed later in Section 3.5.
From the CASSCF data in Tables 3-5, the following conclusions can be drawn: (1) In a close analogy to molecular hydrogen, the matrix elements of polarisability anisotropy are affected by the harmonic approximations to a larger degree than the matrix elements of mean polarisability. This observation concerns all of the studied molecules, with larger deviations detected for lighter molecules. (2) Within a pair of isotopologues (i.e. HF and DF, HCl and DCl, CO and 13 C 16 O, N 2 and 14 N 15 N), the deviations between the matrix element computed using the exact polarisability invariants and the matrix element computed using approximate polarisability invariants are similar, despite that the matrix elements for both for isotopologues may differ to a large degree. This is not surprising since the pairs of isotopologues share the same polarisability invariants. Consequently, similar observation concerns also the total Raman intensities of isotopologues. (3) In the interval [r min , r max ], corresponding to the non-negligible values of the ro-vibrational wave functions, the differences between the exact polarisability invariants (ᾱ (exact) and γ (exact) ) and their linear approximants (ᾱ (1) and γ (1) ) are rather small for all of the studied molecules except F 2 , for which the shape of both the invariants resembles a convex hull in this region. Consequently,ᾱ (1) and particularly γ (1) are bad approximations toᾱ (exact) and γ (exact) , respectively, which results in large deviations between the exact and approximate matrix elements of polarisability invariants and the resulting total Raman intensities. (For a relevant comparison between F 2 and N 2 , see Figure 5.) (4) The deviations from the exact results caused by the double harmonic approximation-corresponding to the choice of ψ harmn andᾱ (1) and γ (1) in the integral in Equation (10)-are molecule dependent and range between −7.5% and +7.5% for theᾱ matrix elements, between −9.5% and +17.5% for the γ matrix elements, and between −8.2% and +9.5% for the total Raman intensities. The deviations are Figure 5. A comparison betweenᾱ exact andᾱ (1) (left panels) and between γ exact and γ (1) (right panels) for F 2 (upper panels) and N 2 (lower panels), together with the relevant vibrational wave functions ψ exact v=0,J=0 and ψ exact v=1,J=0 depicted in the background in grey to illustrate the extent of the integration region in the integral in Equation (10). It is clear that for N 2 , the exact polarisability tensor invariants α exact and γ exact are well modelled by their linear approximantsᾱ (1) and γ (1) , while for F 2 , such an approximation leads to significant errors. The polarisability invariants are computed using CASSCF.
slightly smaller for the heavier isotopologue within a given pair of isotopologues. (5) The accuracy of the single harmonic approximation is 30-40% better than the accuracy of the double harmonic approximation. An exception for this rule is F 2 , for which both approximations are very bad (error of > 17% for γ and > 7% forᾱ). Typical errors attributed to the harmonic approximations are close to several percents.

Uncertainties arising from inconsistencies in the quantum chemical description of polarisability invariants
In Section 3.1.1 we have demonstrated that the distancedependent polarisability tensor components α and α ⊥ computed with various ab initio techniques can differ substantially due to the insufficient description of the electronic wave functions in the vicinity of the dissociation limit by the single-reference-based methods; for details, see Figure 4. The discrepancies between the computed curves are sizable even in the integration interval [r min , r max ] relevant to Equation (10), where the four tested ab initio techniques display considerable disparities (for details, see Figures 2 and 3). It is not surprising that the polarisability matrix elements computed with those curves and tabulated in Tables 3 and 4 differ substantially and affect the resulting Raman intensities tabulated in Table 5. Consequently, large differences in the values of the total Raman intensities computed with different theoretical methods (HF, CASSCF, DFT and CCSD(T)) are observed. A few regularities (and irregularities) are worth highlighting: (1) The polarisability invariants matrix elements computed with the exact methodology with various ab initio techniques may differ substantially. For HF, HCl, CO, and N 2 and their isotopologues, the detected differences can be as large as 56% for the mean polarisability and 60% for the polarisability anisotropy; both these maximal deviations are observed between the matrix elements of N 2 derived from the CASSCF and Hartree-Fock polarisabilities. (2) These maximal deviations have been determined without taking into account the F 2 molecule, for which the deviations are enormous: the matrix elements determined using the Hartree-Fock polarisability invariants are over 15 times larger (ᾱ) and over 60 larger (γ ) than the corresponding CASSCF quantities, while for DFT and CCSD(T) these factors are approximately 6 and 18 times larger, and 3 and 6 times larger, respectively. (3) To understand the strong dependence of the computed matrix elements and the total Raman intensities of F 2 on the selected ab initio technique, we have plotted in Figure 6 a close-up view of the distancedependent polarisability invariantsᾱ exact and γ exact Figure 6. A comparison between distance-dependent polarisability invariantsᾱ exact (left panels) and γ exact (right panels) for F 2 (upper panels) and N 2 (lower panels) computed with various ab initio methods, together with the relevant vibrational wave functions ψ exact v=0,J=0 and ψ exact v=1,J=0 depicted in the background in grey to illustrate the extent of the integration region in the integral in Equation (10). It is clear that for N 2 , all the methods produce curves that can be approximated by linear functions with similar slopes, while for F 2 , distinct functional dependence of curves obtained with different methods causes large differences in the computed matrix elements and total Raman intensities.
computed with various methods for F 2 and N 2 . The difference between these two molecules is striking. For N 2 , all the curves inside the integration interval [r min , r max ] could be reasonably well approximated by linear functions with similar slopes. For F 2 , similar operation could be performed only for the DFT and HF curves, while the curves computed with both CASSCF and CCSD(T) have strong quadratic character. This observations explains well why F 2 behaves differently than other molecules. (4) There seem to be no obvious correlation between the method of calculations, the molecule under study, and the magnitude and the direction of the deviation from the CASSCF results. For every molecule and for every method, the deviations can be either small or large and either positive or negative. (5) Within a pair of isotopologues (i.e. HF and DF, HCl and DCl, CO and 13 C 16 O, N 2 N 2 and 14 N 15 N), the deviations between the matrix elements computed with various ab initio techniques display similar trends. Similar observation concerns also the total Raman intensities of isotopologues.
The deviations from CASSCF in the matrix elements computed with DFT, HF, and CCSD(T) automatically propagate to the associated deviations in the total Raman intensities. These deviations are illustrated graphically in Figure 7, where for each pair of isotopologues, we have represented the CASSCF total Raman intensity as a blue bar of height 100 and the total Raman intensities computed with other ab initio methods as bars in other colours, each of a height associated with the relevant total Raman intensity corresponding to a given method.
The most important observations concerning these results can be summarised as follows: (1) For HF, HCl, CO, and N 2 and their isotopologues, the difference between the reference CASSCF total Raman intensities and analogous results computed with other methods can be quite substantial with deviations from 5% to 57%. For F 2 , the differences are much larger: 264% for CCSD(T), 793% for DFT and over 2000% for HF. (2) The deviations associated with using inadequate electronic structure method to compute distancedependent polarisability invariants completely overshadow the inaccuracies associated with double harmonic approximation. It seems at the moment that solving the problem how to produce accurate distance-dependent polarisability invariants is much more pressing than addressing other methodological details in the process of determination of accurate total Raman intensities.

Discussion
The tested here double harmonic approximation is characterised by errors between −8.2% and +9.5% for the 12 studied molecules. This is not small, because typical experimental errors associated with the acquisition of Raman intensities are lower and could be estimated to be within ±5% [75]. Therefore, precise Raman investigations aiming at accuracy better than 5% should not use double harmonic approximation for the interpretation of the experimental results. An alternative to the standard approach is the scheme presented by us in the current manuscript. However, the actual correspondence between the theoretical off-resonance Raman cross-sections and the experimental Raman intensities is more complicated. The reported here total Raman intensities correspond to the frequency-independent Raman cross-sections, ( 0 /π ) 2ν−1 0ν −3 s × (dσ/d ), expressed in atomic units, withν 0 denoting the absolute frequency of the incident light andν s denoting the absolute frequency of the scattered light [76,77]. A direct comparison of these quantities with experimentally determined Raman crosssections would involve the inclusion of the frequency factorν 0ν 3 s , which takes into account the information about the specific laser used in the experiment and the particular molecular transition. Moreover, experimental Raman intensities include also the effect of temperature affecting the Boltzmann population of the initial state, and the wavelength-dependent response function of the Raman spectrometer. The first of this quantities is relatively straightforward to account for, while the determination of the correction associated with the wavelengthdependent sensitivity of the Raman spectrometer can be by itself an arduous task [78,79].
In a sense, the accuracy estimates of the double harmonic approximation given in the previous paragraph fulfil the main goal of the current study. However, in Figure 8. Comparison of errors originating from the double harmonic approximation (upper panel) and the single harmonic approximation (lower panel) in the process of determination of total Raman intensities for the fundamental vibrational transition for 12 diatomic molecules. The errors correspond to the differences between results computed with the exact methodology based on anharmonic rovibrational wave functions ψ anhrm and distance-dependent polarisability invariantsᾱ exact and γ exact , and the approximate results based on the harmonic ro-vibrational wave functions ψ harmn and linear approximantsᾱ (1) and γ (1) (double harmonic approximation) or on the anharmonic ro-vibrational wave functions ψ anhrm and linear approximantsᾱ (1) and γ (1) (single harmonic approximation).
the process of constructing the answer to the question risen at the beginning of this manuscript, we have realised that there exists much more pressing problem that needs to be urgently addressed in order to provide the physical chemistry community with a method of producing reliable theoretical estimates of total Raman intensities. Namely, we have discovered that the currently available quantum chemical techniques of determination of distance-dependent molecular polarisabilities suffer from serious drawbacks associated with inadequate description of molecular electronic structure at larger internuclear distances. It is well known that single reference methods cannot describe adequately molecules with stretched bonds, which is associated with their inability of a correct description of the dissociation limit. We have discovered here that this limitation affects already the accuracy of molecular polarisability components determined with HF, DFT, and CCSD(T) for molecules with bonds stretched out of the equilibrium position by a relatively small amount. The resulting distancedependent polarisability components computed using HF, DFT, and CCSD(T) differ significantly from each other and diverge systematically from the correct longrange behaviour. We have circumvented this problem by applying in the current study the CASSCF method to determine the polarisability invariant, because it seems at the moment to be the only method able to correctly describe the physically sane behaviour of polarisability invariants at the whole range of the internuclear separations relevant in the process of computing total Raman intensities. However, CASSCF does not account for dynamic correlation and consequently, the polarisability components determined in this way are not too accurate. We hope that this explicit demonstration of the problems associated with determination of accurate distance-dependent polarisability components will stimulate the community for finding better and more reliable methods of computing this quantity. In this regard, CASPT2 (second-order complete active space perturbation theory) and NEVPT2 (second-order n-electronic valence state perturbation theory) techniques [27] would be the most promising way of computing physicallymeaningful distance-dependent polarisabilities, but at the moment the CASPT2 and NEVPT2 response theory codes are not yet available.

Conclusions
We show in this study that the total Raman intensities of fundamental vibrational transitions for 12 diatomic molecules computed within the double harmonic approximation are characterised by errors between −8.2% and +9.5% with respect to an exact treatment based on anharmonic ro-vibrational wave functions ψ anhrm and distance-dependent polarisability invariants α exact and γ exact . A graphical compilation of the resulting errors is given in Figure 8 for all the studied here molecules. These errors are slightly larger than the analogous errors originating from the single harmonic approximation. The largest discrepancies are observed for F 2 , for which the exact functional dependence of the polarisability invariantsᾱ exact and γ exact cannot be well approximated by linear functionsᾱ (1) and γ (1) , respectively.
For isotopologues of molecular hydrogen, the reported here values of total Raman intensities should be close to the exact absolute values. The situation is quite different for heavier molecules, HF, DF, HCl, DCl, CO, 13 C 16 O, N 2 , 14 N 15 N, and F 2 , for which we have realised (unexpectedly) that the determination of accurate distancedependent polarisability invariants needed to evaluate the integral in Equation (2) is difficult if not impossible with the currently existing quantum chemical codes. The spatial extent of the ro-vibrational wave functions requires the knowledge of polarisability invariants at stretched molecular geometries, where the single reference methods start to fail (producing large numerical errors) and the available multireference techniques are not accurate enough (producing likely sizable errors). The errors associated with inaccuracies in the distancedependent polarisability tensor seem to be much larger than the errors originating from the shortcomings of the harmonic approximations, and at the moment seem to be principal factor determining the accuracy of the computed total Raman intensities. The best example here is the F 2 molecule, for which the total Raman intensities computed with various quantum chemical methods can differ even by a factor of 20 (for details, see Section 3.4 and Figure 7). These findings clearly show that there exist an urgent need for developing more accurate methods of computing the distance-dependent polarisability invariants, applicable also to distorted molecular geometries. We believe that this discovery is the most important outcome of our work. Notes 1. Total wave function here is expressed as product of harmonic oscillator wave functions of each normal mode, i.e. ψ v i = φ v i (Q k ), and similarly ψ v f = φ v f (Q k ). 2. This approach also covers the infrared intensities where the polarisability operator is replaced by the dipole moment operator. 3. We limit our discussion to the transitions in the ground electronic state which is non-degenerate. Only offresonance Raman scattering intensities are studied here.