Complete analytic anharmonic hyper-Raman scattering spectra

We present the first computational treatment of the complete second-order vibrational perturbation theory applied to hyper-Raman scattering spectroscopy. The required molecular properties are calculated in a fully analytic manner using a recently developed program [Ringholm, Jonsson and Ruud, J. Comp. Chem. , 2014, 35 , 622] that utilizes recursive routines. For some of the properties, these calculations are the first analytic calculations of their kind at their respective levels of theory. We apply this approach to the calculation of the hyper-Raman spectra of methane, ethane and ethylene and compare these to available experimental data. We show that the anharmonic corrections have a larger eﬀect on the vibrational frequencies than on the spectral intensities, but that the inclusion of combination and overtone bands in the anharmonic treatment can improve the agreement with the experimental data, although the quality of available experimental data limits a detailed comparison.


Introduction
Hyper-Raman scattering was predicted in the late 1950s 1 and observed experimentally in the mid-1960s. 2 It is a process which complements the infrared (IR) and conventional Raman spectroscopies because the vibrational transition moments involve the molecular first hyperpolarizability instead of the dipole moment or polarizability tensors that govern IR and Raman spectroscopy, respectively. Vibrational transitions that do not give an IR or conventional Raman signal could therefore give a hyper-Raman signal. [3][4][5] Conversely, congested spectral features in IR or Raman spectroscopy could become less congested in hyper-Raman spectroscopy. Moreover, hyper-Raman spectroscopy has been found to be particularly useful experimentally for probing the frequency region associated with low-energy vibrational transitions.
Spectroscopies involving molecular vibrations have been simulated ab initio for a long time, and the stage was set by early developments that allowed for the calculation of geometric derivatives of the molecular energy. [6][7][8][9][10] Subsequent developments have made the calculation of second-order geometric derivatives available for various levels of electronic structure theory. 8,[10][11][12][13][14][15][16][17][18][19][20][21] The most commonly used model for simulating vibrational spectroscopies is the so-called double-harmonic approximation, 22,23 where the vibrational potential between atoms is described by a harmonic oscillator model, and where the spectroscopic intensities are governed by the first geometric derivatives of the relevant polarization properties. For hyper-Raman scattering, this polarization property is the molecular first hyperpolarizability, and calculations of this property and its subsequent use in theoretical hyper-Raman spectra has been carried out previously, in particular by the pioneering work of Quinet, Champagne and co-workers, [24][25][26][27][28] but also more recently by groups involving some of the present authors. 5,29 The double-harmonic approach, however, ignores anharmonicity in the vibrational motions and disregards spectral effects stemming from excitations involving more than one vibrational quantum, and the inclusion of anharmonic effects will therefore tend to improve the quality of the computed results. [30][31][32][33][34] Including anharmonic effects in spectroscopies involving vibrations is then a twofold task: correcting the vibrational frequencies and the intensities. Correcting vibrational frequencies by perturbation theory has been done for quite some time, 32,[35][36][37][38] whereas the development of corrections for the intensities is more recent and has been extensively explored by e.g. Barone and Bloino for IR and Raman spectroscopy. 34 Except for a partial treatment in an earlier study, 26 the calculation of anharmonic corrections in hyper-Raman scattering has not been reported. In this work, we present the first complete second-order vibrational perturbation theory (VPT2) treatment of hyper-Raman scattering.
The properties that are needed in order to evaluate the expressions that the VPT2 treatment generates (for both frequencies and intensities) are of higher order and their calculation is challenging. In particular, for hyper-Raman scattering, a full second-order treatment requires the calculation of properties such as the third-and fourth-order geometric derivative of the energy, henceforth called the cubic and quartic force constants, respectively, and the third-order geometric derivative of the molecular first hyperpolarizability. The latter property can be formulated as a sixth-order derivative of the molecular quasienergy and its calculation, because of its high order, is complicated due to the numerous contributions that are involved, as discussed previously for the case of fifth-order derivatives. 39 The calculation of the properties needed for anharmonic corrections can be done by numerical or analytic methods through the use of response theory, but as is well known 40 and as demonstrated in recent work 39 for IR and Raman spectroscopy, calculations involving numerical methods can be prone to errors and it can be difficult to identify when these errors actually have occurred. For these reasons, we recommend and prefer analytic methods whenever available. However, analytic calculations have traditionally been considered to be substantially more challenging to implement than numerical methods, and the majority of studies to date have used numerical methods for these purposes, although there exist some computer programs that permit calculations of some of the necessary properties. 41,42 We have recently created an implementation 43 of an openended formulation of response theory 44 that can manage the analytic calculation of any response property, provided that external program routines [45][46][47] can supply various contributions that may then be assembled into the desired molecular property. The implementation uses recursive routines to ensure generality and open-endedness, and we have already used it for the first analytic calculations of many properties of relevance to vibrational spectroscopies at the density-functional theory (DFT) level. This includes the first reported analytic DFT quartic force constants 48 and geometric first derivatives of the hyperpolarizability for hyper-Raman scattering in the doubleharmonic approximation. 29 We have also recently applied it for a complete analytic DFT VPT2 treatment of IR and Raman spectroscopy, 39 and in the present work, we extend this to hyper-Raman spectroscopy. In all our previous work as well as in the present work, calculating the properties needed, including the sixth-order property mentioned above, amounts to a routine execution of our program, and in light of the problems that may be introduced by the use of numerical methods, we believe that this constitutes a significant improvement over previous approaches. We note, however, that due to current program limitations, the analytic calculation of the sixth-order derivatives are limited to the Hartree-Fock level of theory. All other derivatives up to and including fifth-order energy derivatives have been calculated analytically using DFT.
The rest of this paper is organized as follows: in Section 2, we present the theory of hyper-Raman scattering and the response theory underlying the calculation of the requisite response properties. In Section 3, we show the computational details of our calculations, the results of which are presented and discussed in Section 4. Finally, we give some concluding remarks and an outlook in Section 5.

Theory
For completeness, we give in Section 2.1 an introduction to the theory of hyper-Raman scattering that follows closely the presentation in our previous work. 29 In Section 2.2, we show how our approach can be used to calculate anharmonic corrections to hyper-Raman spectra. Finally, we present briefly in Section 2.3 how our open-ended implementation of response theory can be used to calculate the properties that are needed to evaluate the expressions, including VPT2 anharmonic corrections, that determine the spectroscopic intensities and vibrational frequencies.

Hyper-Raman
A sample exposed to an incident laser field of intensity I 0 and frequency o will produce, for a vibrational transition of frequency o a , a Stokes hyper-Raman band at frequency (2o À o a ) of intensity 24,26,28 where the initial and final states of the system, with population N, have been denoted |C init i and |C fin i, respectively. The first hyperpolarizability is denoted as b(o, Q) at frequency o and is a function of the nuclear displacements along the reduced (frequency-weighted) normal coordinates Q.
If the Born-Oppenheimer approximation is applied, the different states can be written as a product of pure vibrational and electronic states. However, this neglects the vibronic coupling which may be important close to an electronic resonance. For (2o À o a ) detuned from electronic resonance, the electronic state of both the initial and final states can be taken to be the electronic ground state. It is then possible to integrate over the electronic degrees of freedom, giving the intensity as with |ni and |pi being, respectively, the initial and final vibrational states and b the electronic first hyperpolarizability of the ground state.
When considering an incident plane wave linearly polarized in the a direction, the total non-polarized scattered hyper-Raman intensity can be further simplified and is then given by: 49 where k denotes the Boltzmann constant and T the temperature, entering in the Boltzmann term in order to take into where b 0 abg is used here to denote hb abg i 0i . These expressions are obtained by rearranging eqn (10) and (11) in ref. 49.
The hyper-Raman intensity can then easily be split into the contribution from the vertically (VV) and horizontally (HV) polarized scattered light:

Anharmonic corrections to vibrational frequencies and spectroscopic intensities
The first hyperpolarizability and the vibrational potential can be expanded around the equilibrium geometry along the reduced normal coordinates Q i of the system in a Taylor series, giving where Greek subscripts denote Cartesian coordinates, and Note that, since the expansion is carried out at the equilibrium geometry, qV eq /qQ i | eq = 0 for every reduced normal coordinate i.
In order to ease the reading of the equations, the compact notation introduced on the last right-hand side of each of eqn (8) and (9) will be used in the rest of this section. The most commonly used approach in the study of spectroscopies involving molecular vibrations is the so-called doubleharmonic approximation. In this approximation, the molecular vibrations are described by harmonic oscillators and the spectroscopic intensities are calculated from the first geometric derivative of the relevant polarization properties, and the expansions in eqn (8) and (9) are therefore truncated at first and second order, respectively. The transition moment of the first hyperpolarizability, to be used in eqn (3)-(7), for a transition from the ground state to the vibrationally excited state corresponding to the ith normal mode, will then be given by b abg Once the harmonic vibrational frequencies and normal modes of vibration have been determined from the well-established eigenanalysis of the molecular Hessian, 23 a second-order perturbational approach can be used to obtain anharmonic corrections to the fundamental vibrational frequencies and frequencies corresponding to the two-quantum combination and overtone excited vibrational states. When using Vibrational Perturbation Theory to second order (VPT2), 30-32 the harmonic singly excited normal mode frequencies o i , the cubic V (3) ijk and quartic V (4) ijkl force constants, the rotational constants B a and the Coriolis vibrationrotation coupling constants z a ij are used to compute the corrected anharmonic vibrational frequencies of the normal modes n i 1, of the first overtones n i 2 and the first combination bands n i 1 j 1.
n i 1 j 1 = n i 1 + n j 1 + X ij (13) where the terms X ii and X ij are given as and with O ijk defined as In this work, the Generalized Vibrational second-order Perturbation Theory (GVTP2) model 34,38,50 will be used, in which the standard VPT2 scheme is used for the molecular vibrations when there is no Fermi resonance, but where the terms in the VPT2 treatment that are affected by Fermi resonances are not included 32 and are instead treated variationally.
The spectral intensities can also be corrected by a similar perturbational approach. The GVPT2 expressions for IR and Raman intensities have already been derived and presented in earlier works by Bloino and Barone. 34,50 In this work, the same approach is applied to hyper-Raman intensities, and the transition moment of the first hyperpolarizability in the GVPT2 approach can then be obtained from the following expressions for a transition from the ground state to the vibrationally excited state corresponding to the ith normal mode (hb abg i 0,i 1), to the first overtone of the ith normal mode (hb abg i 0,i 2) or to the combination of the ith and the jth normal modes (hb abg i 0,i 1 j 1):

Analytic calculation of response properties
This section will introduce the response theory that underlies the recursive implementation from which the high-order molecular properties that enter into the expressions presented in Section 2.1 are calculated. This method has been presented thoroughly in the original work 43,44 and we will therefore only give a brief outline here.
A response property corresponding to a collection of perturbations a, b, c. . . can be expressed as derivatives of a generalized quasienergy Lagrangian L abc k,n , where the superscript denotes perturbation-strength differentiation and subsequent evaluation at zero perturbation strength. The general expression for a response property, with collected associated frequencies of the perturbations b, c. . . defined as o bcÁ Á Á , is given by where ¼ fTrg T symbolizes that the trace and time average of the right-hand side is taken, and where the choice of integers (k, n) and the various forms in which they occur in the subscript determines the truncation order for the perturbed arguments, notably the perturbed or unperturbed density matrix D or generalized Kohn-Sham Fock matrixF that enter into the terms on the right-hand side, so that a generalization is made from the well-known (n + 1) and (2n + 1) rules to a (k + n + 1) framework, 51 permitting rule choices that are intermediate between the (n + 1) and (2n + 1) extremes. Representing the atomic orbital (AO) overlap matrix as S, and using a tilde to b abg represent a quantity considered at an arbitrary perturbation strength, the generalized Fock matrix and energy are given bỹ and respectively, where the constituent terms are the one-electron matrix h and two-electron matrix G g with fractional exchange g A [0,1], the external field operator Ṽ t , the half-timedifferentiated overlap matrix T, the exchange-correlation contributions F xc and Ẽ xc [r(D)] for a density r, and the nuclearnuclear repulsion term h nuc . Also introduced in eqn (20) are the energy-and frequency-weighted Fock matrix W , defined as the idempotency condition matrix Ỹ as the time-dependent self-consistent field (TD-SCF) condition matrix Z asZ and their respective Lagrange multiplier termsl a andz a as andz a ¼F aDS À where we have introduced the notation and where adjungation is defined to take place before any time differentiation. Eqn (20) can be evaluated and the various perturbed density and Fock matrices can be determined by making use of a set of recursive algorithms, and we refer to the original work 43,44 for details about this procedure. The benefits of a recursive approach can be demonstrated by considering the highest-order property calculated in this work, namely the third-order geometrical derivative of the first hyperpolarizability.
We note in passing that in order to evaluate this expression with the (k, n) = (2, 3) rule choice, it is necessary to calculate perturbed Fock and density matrices up to the combined second order of geometrical perturbation and first order of electric dipole perturbation. A part of the procedure to obtain these matrices, and typically the most computationally demanding step, is solving a set of linear equation systems for the so-called response parameters, 44 and, taking M to be the number of geometrical coordinates of the system under study, the number of such equation systems to be solved is proportional to M 2 . The term is given by where we have written terms that collectively are permutations of identical operators only once and removed terms that vanish upon differentiation because of lack of dependence on the perturbation operators, and where we introduced the notation and where tracing is exemplified by Using the same convention as was done for eqn (31), the other terms in eqn (30) can be written as and The complexity of these expressions, in particular eqn (31), coupled with the fact that constituent terms in eqn (34)-(37) such as W ggf À2o fofo 3 0 may each contain a large number of contributions, would make the creation of a tailored program routine for the calculation of this property a very cumbersome task, justifying the use of a recursive approach as done in this work. We remark that the exchange-correlation contribution E gggf À2o fofo xc in eqn (31) has not yet been implemented in the presently used program module, 52 and we have for this particular property therefore performed the calculations at the HF level of theory in the manner described in Section 3.

Computational details
The various tensors required for the calculation of the vibrational spectra, with the exception of the Hessian, have been computed using the recursive implementation 43 of the open-ended response theory framework of Thorvaldsen et al., 44 which has been implemented in a local version of the DALTON2013 program package. 53,54 Response equations have been solved using the linear response solver of Jørgensen et al. 55 Differentiated oneand two-electron integrals were computed using the GEN1INT 45,56 and CGTO-DIFF-ERI 46,57 programs, respectively, except for some of the lower-order two-electron geometric derivatives which were computed using existing functionality in DALTON. The differentiated exchange-correlation (XC) energy and potential contributions needed in the DFT calculations were computed using the XCFUN library, 47,58 where the integrator XCINT 52 has been used for the integration of the XC contributions. Coriolis coupling constants for anharmonic corrections to the vibrational frequencies are not calculated in a response theory framework, but in the manner outlined in ref. 59.
The calculations have been performed at the DFT level of theory using the B3LYP hybrid functional, 60-62 except the third geometry derivative of the first hyperpolarizability E gggf À2o fofo xc À Á which has been calculated at the Hartree-Fock level due to technical limitations in XCINT. The B3LYP functional has already been shown to give good results for the calculation of higher-order properties in earlier work. 29, 48 Dunning's augmented correlation-consistent polarized triple-z (aug-cc-pVTZ) basis set 63 has been used. The study was conducted for methane (CH 4 ), ethane (C 2 H 6 ) and ethylene (C 2 H 4 ). All calculations have been carried out without any use of point group symmetry in the calculations. The higher-order molecular properties have been computed at the optimized geometry using the recursive response implementation. 43 The molecular Hessian was used in a vibrational analysis to find the harmonic vibrational frequencies and to transform the geometrically differentiated property tensors from a Cartesian basis to a reduced normal coordinate basis 64 to calculate anharmonically corrected frequencies and spectral intensities. We remark that this DFT/B3LYP Hessian was also used for the transformations of the third geometric derivatives of the first hyperpolarizability, for which the Cartesian basis calculations were carried out at the HF level.
Anharmonic corrections to the fundamental frequencies, as well as first overtone and combination band frequencies, were calculated from the cubic and quartic force constants, the rotational constants and the Coriolis coupling constants using the aforementioned GVPT2 model. 34 The spectral band shapes have been modelled using Lorentzian functions with a 10 cm À1 width at half maximum. A 1 cm À1 resolution was used to plot all spectra.
First-order geometric derivatives of the respective polarization properties in reduced normal coordinate basis were used for the evaluation of the harmonic spectral intensities for the different vibrational spectroscopies. Anharmonic corrections to the intensities were then calculated using also the second-and third-order geometric derivatives of the corresponding properties and the cubic and quartic force constants in a reduced normal coordinate basis, using the GVPT2 model. The threshold for identifying Fermi resonances was set to the default value of 1 cm À1 for the Martin parameter 32 as used by Bloino and Barone. 34

Results and discussion
In this section, the infrared (IR), Raman and hyper-Raman spectra of three structurally related molecules for which experimental hyper-Raman spectra exists have been calculated using the analytic approach. For these test molecules, the experimental hyper-Raman spectra from Verdieck et al. 65

Methane
The calculated harmonic and anharmonic vibrational frequencies for methane are collected in Table 1 together with experimental values. 66 While the methane molecule is a spherical top, the vibrational perturbation theory used in this work is formally valid only for asymmetric top molecules, and the effect of this has been neglected since we believe that it will not have a substantial impact on the results. Because of the T d symmetry of methane, several vibrational modes are degenerate and thus only 4 distinct vibrational frequencies corresponding to fundamental excitations are observed experimentally. These degeneracies are well reproduced by our theoretical calculations at the harmonic level with one mode at 3026 cm À1 , two at 1556 cm À1 , three at 3126 cm À1 and three at 1339 cm À1 . These degeneracies are conserved when adding the anharmonic corrections, demonstrating the numerical stability of the analytic approach, where symmetry has not been imposed during the calculations. Theoretical harmonic frequencies are in qualitative agreement with the experimental fundamental frequencies.
The agreement is improved when applying anharmonic corrections. The improvement goes from a deviation larger than 100 cm À1 to less than 30 cm À1 for the two high-frequency modes and from a deviation larger than 20 cm À1 to 13 and 1 cm À1 for the two lowfrequency modes, respectively. Fig. 1 shows the calculated infrared and non-polarized Raman (Fig. 1a) and non-polarized hyper-Raman (Fig. 1b) spectra of methane at the harmonic and anharmonic level. The experimental hyper-Raman spectrum from ref. 65 is also shown for comparison in Fig. 1c.
For symmetry reasons, two peaks are expected to be observed in IR, four in Raman and three in hyper-Raman in the double-harmonic approximation, as also observed in our theoretical harmonic spectra. In the case of the Raman spectrum, the peak 4 1 that is expected at 1339 cm À1 is predicted, but with an intensity so low compared to the others that it does not appear in Fig. 1a.
In both the harmonic and anharmonic theoretical hyper-Raman spectra, three peaks are present, corresponding to modes 1 1 , 3 1 and 4 1 . The 3 1 peak is the most intense in both cases. The main improvement going to the anharmonic treatment lies in the correction of the vibrational frequencies, as was already shown in our previous work for the IR and Raman spectra. 39 However, including anharmonic corrections increases the difference in intensity between the strong 3 1 peak and the two other weak ones. In the experimental hyper-Raman spectrum, one feature appears larger than the experimental noise around 3000 cm À1 . This feature corresponds to the intense 3 1 peak in the theoretical spectra. A shoulder on the lower-frequency side of the feature (around 2900 cm À1 ) appears in the experimental spectrum, and can be linked to the 1 1 peak predicted in the theoretical spectra. In the experimental reference article, 65 a peak is also mentioned at 1306 cm À1 , corresponding to 4 1 . As can be seen in Fig. 1c, this peak is so weak that it is almost impossible to distinguish it from the experimental noise. This low-intensity peak is well reproduced in both the harmonic and anharmonic theoretical spectra. Due to the difficulty of distinguishing the 1 1 and 4 1 peaks from the noise in the experimental spectrum, it is not possible to conclude whether the harmonic or anharmonic treatment better reproduces the experimental relative intensities, only that our theoretical treatment reproduces well the largest intensity of the 3 1 peak.

Ethane
The calculated harmonic and anharmonic vibrational frequencies for ethane are compiled in Table 2 and compared to experimental values. 66 Because of the D 3d symmetry of ethane, several vibrational modes are degenerate and thus only 12 distinct vibrational frequencies corresponding to fundamental excitations are observed experimentally. Modes 7,8,9,10,11 and 12 are all doubly degenerate. These degeneracies are well reproduced by our theoretical harmonic calculated results and are conserved when adding the anharmonic corrections. However, the harmonic treatment gives an additional pseudodegeneracy between modes 1 1 and 5 1 (3024 cm À1 ), both of these modes corresponding to symmetric stretching motion of the CH 3 moieties (one fully symmetric, the other antisymmetric with respect to the inversion). The anharmonic treatment removes this accidental degeneracy (1 1 at 2868 cm À1 and 5 1 at 2877 cm À1 ). While experimentally 1 1 is higher in frequency than 5 1 , the anharmonic treatment predicts the opposite. Globally, the anharmonic frequencies are in better agreement with the experimental spectra, with the notable exception of modes 1 1 , 3 1 and 12 1 with absolute deviations from the experimental value going from 70 to 86 cm À1 , from 1 to 24 cm À1 and from 3 to 4 cm À1 , respectively. Apart from these three modes, the harmonic treatment shows an overestimation of the vibrational frequencies compared to experiment of more than 90 cm À1 for the modes above 2000 cm À1 , between 30 and 35 cm À1 for the modes in the 1000 to 2000 cm À1 frequency range, and less than 15 cm À1 for the lower-frequency modes. Adding anharmonic corrections reduces the absolute deviation from experiment to less than 40 cm À1 for the modes above 2000 cm À1 and less than 10 cm À1 for the modes vibrating at lower frequencies. Fig. 2 shows the calculated infrared and non-polarized Raman (Fig. 2a) and non-polarized hyper-Raman (Fig. 2b) spectra of ethane at the harmonic and anharmonic levels. The experimental hyper-Raman spectrum from ref. 65 is also reproduced in Fig. 2c for comparison.
For symmetry reasons, five peaks are expected in IR, six in Raman and six in hyper-Raman in the double-harmonic approximation. This partitioning of the peaks is well reproduced in our theoretical results. In the hyper-Raman spectrum, all six peaks are visible, but in the Raman and IR spectra, some peaks (6 1 at 1411 cm À1 for IR and 2 1 at 1420 cm À1 and 9 1 at 1221 cm À1 for Raman) are predicted with intensities so low compared to the highest-intensity peak that they are not visible in the spectra presented here. The 4 1 peak at 302 cm À1 is present in the hyper-Raman spectrum but neither in the IR nor in the Raman ones.
The most notable improvement brought by the addition of anharmonic corrections in these three spectra lies again in the correction of the frequencies. Relative spectral intensities are not strikingly modified by the addition of anharmonic corrections. The main modification in terms of the intensity in the hyper-Raman spectrum is the increase of the intensity of the highest-frequency peak corresponding to the 10 1 mode. Another interesting improvement brought by the anharmonic treatment is the inclusion of peaks corresponding to combination modes. The most intense combination peak is found at 2914 cm À1 , corresponding to the 8 1 11 1 combination band. Two very weak peaks can also be seen at 2745 and 3755 cm À1 corresponding, respectively, to the 2 1 6 1 and 7 1 12 1 combination Table 2 Theoretical harmonic o i , anharmonic n theo i and experimental n expt i normal mode vibrational frequencies of ethane (in cm À1 ). Activity in infrared (IR), Raman (R) and hyper-Raman (hR) spectroscopies given for D 3d symmetry. Symmetry character G and degeneracy n deg are also given for each mode  modes. In the experimental spectrum, four features can be distinguished: a strong band around 3000 cm À1 and three weak ones that are hardly discernible from the experimental noise around 1400, 800 and 300 cm À1 . The main feature corresponds to the group of peaks predicted by the theoretical approach in this region (5 1 , 10 1 and 8 1 11 1 ). The three other peaks correspond respectively to the two-peak group 6 1 and 11 1 , to the 12 1 peak and to the 4 1 peak. The global intensity of the group of peaks around 3000 cm À1 is higher than the intensities of the three other peaks in a similar manner as the main feature at the same position in the experimental spectrum is more intense than the other features. In both the harmonic and anharmonic calculated spectra, the global intensities of the three lowerintensity features are similar to each other.

Ethylene
The calculated harmonic and anharmonic vibrational frequencies for ethylene are compiled in Table 3 and compared to experimental values. 66 Because of the D 2d symmetry of ethylene, five peaks are expected in the IR spectrum, six in the Raman, and six in the hyper-Raman spectrum at the double-harmonic level of approximation. This partitioning of the peaks is again well reproduced by our theoretical approach. In the calculated Raman spectrum, only five peaks seem to appear, but a sixth one, the 6 1 peak at 1245 cm À1 , is also predicted but with a relative intensity so low that it does not appear in the spectrum.
In hyper-Raman, the 12 1 peak at 1478 cm À1 also has a very low intensity and is almost not visible in the theoretical spectrum. The 4 1 peak at 1060 cm À1 is present in the hyper-Raman spectrum, but neither in the IR nor in the Raman spectra. The harmonic treatment overestimates vibrational frequencies by more than 30 cm À1 for all modes but 6 1 and 10 1 , which are overestimated by less than 10 cm À1 . For the modes with the higher frequencies (above 2900 cm À1 ), the overestimation is larger than 85 cm À1 . All anharmonic frequencies are in agreement with the experimental ones, with absolute deviations lower than 30 cm À1 , with the exception of mode 5 1 for which the deviation is 51 cm À1 . Only the two normal modes with harmonic deviation less than 10 cm À1 (6 1 and 10 1 ) present a better agreement with experiment at the harmonic level than at the anharmonic level for the vibrational frequencies. However, the absolute deviation worsens only by 1 cm À1 , still remaining lower than 10 cm À1 . Fig. 3 shows the calculated infrared and non-polarized Raman (Fig. 3a) and non-polarized hyper-Raman (Fig. 3b) spectra of ethylene at the harmonic and anharmonic levels. The experimental hyper-Raman spectrum from ref. 65 is reproduced in Fig. 3c for comparison. Once again, the main improvement from the anharmonic treatment compared to the harmonic one comes from the frequencies. Relative intensities of low-intensity modes are slightly decreased and relative intensities of high-intensity modes are slightly increased when including anharmonic corrections. Another improvement in the anharmonic spectrum is the appearance of a very weak peak at 2019 cm À1 arising from the 4 1 7 1 combination mode. In the experimental spectrum, two main features can be observed around 1000 cm À1 and around 3000 cm À1 . These two features can be linked respectively to the 4 1 , 7 1 and 10 1 three-peak group and the 9 1 and 11 1 twopeak group. The intensity of the 1000 cm À1 group is higher than the one of the 3000 cm À1 group in a similar way that the 1000 cm À1 feature is more intense than the 3000 cm À1 feature. The two other very weak peaks found in the theoretical spectra are too weak to be distinguished from the noise in the experimental spectrum.
The calculated (harmonic and anharmonic) hyper-Raman spectra using the HV polarization setup is show in Fig. 4a. The corresponding experimental spectrum from ref. 65 is reproduced in Fig. 4b.
In the HV-polarization case, the anharmonic treatment brings the same corrections to the intensities as for the non-polarized hyper-Raman case. In the anharmonic HV-polarized spectrum, a few more weak peaks become visible, such as the 5 1 6 1 combination mode at 2608 cm À1 and the 5 1 7 1 combination mode at 3962 cm À1 . Also, the relative intensities of the two groups of peaks around 1000 cm À1 and 3000 cm À1 in the theoretical spectra are in agreement with the experimental ones, but the dominant contributing mode to the feature around 1000 cm À1 changes from 7 1 to 4 1 when going from the non-polarized to the HV polarization setup, which is not apparent from the experimental results, and this supports the notion that the choice of polarization setups offers another degree of freedom for probing different vibrational modes of the system under study. Finally, we remark that as expected, the quality of our calculated results seems unaffected by the choice of polarization setup.
Overall, we see that the inclusion of anharmonic effects produces an improvement in the quality of the calculated spectra and that this improvement mainly concerns the vibrational frequencies, even though there for these frequencies still exists some disagreement. A significant part of this disagreement is believed to be a consequence of the level of theory chosen for the calculation of the harmonic vibrational frequencies. This has been addressed in earlier works by Barone and co-workers, 67,68  where they establish a computational protocol wherein the larger harmonic contributions to the vibrational frequencies are calculated at a high level of theory, but where the smaller anharmonic corrections are done at a less computationally demanding level. By this approach, the cost of the various contributions are balanced according to their expected magnitude. Analogously, it may be speculated that for the highestorder properties used for anharmonic corrections, calculations at the HF level may offer enough accuracy since the magnitude of the corrections resulting from these properties are expected to be small.

Summary and concluding remarks
We have used our recently developed recursive implementation of response theory to carry out the first fully analytic complete VPT2 treatment of hyper-Raman spectroscopy, including anharmonic corrections to both the vibrational frequencies and the spectroscopic intensities. Using this approach, we have produced IR, Raman and hyper-Raman spectra for methane, ethane and ethylene. The calculations needed to achieve this includes the first analytic ab initio calculations of the respectively second and third geometric derivatives of the first hyperpolarizability at respectively the DFT and HF levels of theory. These properties are respectively fifth-and sixth-order derivatives of the molecular energy, and the task of calculating these and the other properties used in this work is greatly simplified by the use of our recursive approach. In line with previous findings for IR and Raman spectroscopy, 34,39 we observe that the inclusion of anharmonicity has a larger effect on the frequencies than on the intensities, but we  View Article Online also see that the inclusion of combination and overtone bands using the anharmonic treatment can add important features that are not present at the double-harmonic level, such as the 8 1 11 1 combination peak at 2914 cm À1 for the ethane hyper-Raman spectrum. The anharmonic spectra are in qualitatively better agreement with experimental spectra, but for the systems in this work, a more detailed analysis of the improvements obtained by the inclusion of anharmonic effects is hindered by instrumental limitations leading to a significant level of noise in the experimental results. Future development of the analytic approach for the calculation of spectroscopy involving molecular vibrations will aim to include molecular environment and solvent effects through the use of combined quantum-mechanical and classical approaches, and work on this task has already started. 69