Complete analytic anharmonic hyper-Raman scattering spectra

Abstract

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 effect 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.

---

1 Introduction

Hyper-Raman scattering was predicted in the late 1950s¹ and observed experimentally in the mid-1960s.² 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–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–10 Subsequent developments have made the calculation of second-order geometric derivatives available for various levels of electronic structure theory.^8,10–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–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–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–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.³⁴ Except for a partial treatment in an earlier study,²⁶ 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.³⁹

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⁴⁰ and as demonstrated in recent work³⁹ 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⁴³ of an open-ended formulation of response theory⁴⁴ that can manage the analytic calculation of any response property, provided that external program routines^45–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⁴⁸ and geometric first derivatives of the hyperpolarizability for hyper-Raman scattering in the double-harmonic approximation.²⁹ We have also recently applied it for a complete analytic DFT VPT2 treatment of IR and Raman spectroscopy,³⁹ 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.

2 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.²⁹ 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.

2.1 Hyper-Raman

A sample exposed to an incident laser field of intensity [scr I, script letter I]⁰ and frequency ω will produce, for a vibrational transition of frequency ω_a, a Stokes hyper-Raman band at frequency (2ω − ω_a) of intensity^24,26,28

[scr I, script letter I]_hR,a(2ω − ω_a) ∝ Nℏ[scr I, script letter I]²₀(2ω − ω_a)⁴| 〈Ψ_init|β(ω, Q)|Ψ_fin〉|², (1)

where the initial and final states of the system, with population N, have been denoted |Ψ_init〉 and |Ψ_fin〉, respectively. The first hyperpolarizability is denoted as β(ω, Q) at frequency ω 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 (2ω − ω_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

[scr I, script letter I]_hR,a(2ω − ω_a) ∝ Nℏ[scr I, script letter I]²₀(2ω − ω_a)⁴|〈n|[small beta, Greek, macron](ω, Q)|p〉|², (2)

with |n〉 and |p〉 being, respectively, the initial and final vibrational states and [small beta, Greek, macron] the electronic first hyperpolarizability of the ground state.

When considering an incident plane wave linearly polarized in the α direction, the total non-polarized scattered hyper-Raman intensity can be further simplified and is then given by:⁴⁹

where k denotes the Boltzmann constant and T the temperature, entering in the Boltzmann term in order to take into account the population of the excited vibrational states. The averages and can be expressed as:andwhere is used here to denote 〈β_αβγ〉_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:

2.2 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, givingwhere Greek subscripts denote Cartesian coordinates, andNote that, since the expansion is carried out at the equilibrium geometry, ∂V_eq/∂Q_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 double-harmonic 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

Once the harmonic vibrational frequencies and normal modes of vibration have been determined from the well-established eigenanalysis of the molecular Hessian,²³ 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 ω_i, the cubic V⁽³⁾_ijk and quartic V⁽⁴⁾_ijkl force constants, the rotational constants B_α and the Coriolis vibration-rotation coupling constants ζ^α_ij are used to compute the corrected anharmonic vibrational frequencies of the normal modes ν_i¹, of the first overtones ν_i² and the first combination bands ν_i¹j¹.

ν_i² = 2ν_i¹ + 2X_ii, (12)

ν_i¹j¹ = ν_i¹ + ν_j¹ + X_ij (13)

where the terms X_ii and X_ij are given asandwith Ω_ijk defined as

Ω_ijk = (ω_i + ω_j + ω_k)(−ω_i + ω_j + ω_k)(ω_i − ω_j + ω_k)(ω_i + ω_j − ω_k). (16)

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³² 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 (〈β_αβγ〉_0,i¹), to the first overtone of the ith normal mode (〈β_αβγ〉_0,i²) or to the combination of the ith and the jth normal modes (〈β_αβγ〉_0,i¹j¹):

2.3 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 [script 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 ω_bc⋯, is given by

where 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 matrix 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,⁵¹ 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 represent a quantity considered at an arbitrary perturbation strength, the generalized Fock matrix and energy are given byandrespectively, where the constituent terms are the one-electron matrix h̃ and two-electron matrix G̃^γ with fractional exchange γ ∈ [0,1], the external field operator Ṽ^t, the half-time-differentiated overlap matrix T̃, the exchange–correlation contributions F̃_xc and Ẽ_xc[[small rho, Greek, tilde](D̃)] for a density ρ, and the nuclear-nuclear repulsion term h̃_nuc. Also introduced in eqn (20) are the energy- and frequency-weighted Fock matrix W̃, defined asthe idempotency condition matrix Ỹ as

Ỹ = D̃S̃D̃ − D̃, (24)

the time-dependent self-consistent field (TD-SCF) condition matrix Z̃ asand their respective Lagrange multiplier terms [small lambda, Greek, tilde]_a and [small zeta, Greek, tilde]_a as

[small lambda, Greek, tilde]_a = [D̃^aS̃D̃]^⊖, (26)

andwhere we have introduced the notation

[M]^⊖ = M − M^†, (28)

and

[M]^⊕ = M + M^†, (29)

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. Denoting a geometric perturbation by g and an electric dipole perturbation associated with a frequency ω as f_ω, this property can be represented as a perturbation tuple (g, g, g, f_−2ω, f_ω, f_ω). Making the rule choice (k, n) = (2, 3), eqn (20) takes the formWe 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,⁴⁴ 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². The term is given bywhere 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 notationand where tracing is exemplified byUsing the same convention as was done for eqn (31), the other terms in eqn (30) can be written asandThe complexity of these expressions, in particular eqn (31), coupled with the fact that constituent terms in eqn (34)–(37) such as 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 in eqn (31) has not yet been implemented in the presently used program module,⁵² and we have for this particular property therefore performed the calculations at the HF level of theory in the manner described in Section 3.

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⁴³ of the open-ended response theory framework of Thorvaldsen et al.,⁴⁴ 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.⁵⁵ Differentiated one- and 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⁵² 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 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-ζ (aug-cc-pVTZ) basis set⁶³ has been used. The study was conducted for methane (CH₄), ethane (C₂H₆) and ethylene (C₂H₄). 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.⁴³ 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⁶⁴ 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.³⁴ The spectral band shapes have been modelled using Lorentzian functions with a 10 cm⁻¹ width at half maximum. A 1 cm⁻¹ 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⁻¹ for the Martin parameter³² as used by Bloino and Barone.³⁴

4 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.⁶⁵ are shown for comparison and validation purposes. The results for methane, ethane and ethylene are presented and discussed in Sections 4.1, 4.2 and 4.3, respectively.

4.1 Methane

The calculated harmonic and anharmonic vibrational frequencies for methane are collected in Table 1 together with experimental values.⁶⁶ 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⁻¹, two at 1556 cm⁻¹, three at 3126 cm⁻¹ and three at 1339 cm⁻¹. 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⁻¹ to less than 30 cm⁻¹ for the two high-frequency modes and from a deviation larger than 20 cm⁻¹ to 13 and 1 cm⁻¹ for the two low-frequency modes, respectively.

Table 1 Theoretical harmonic ω_i, anharmonic ν^theo_i and experimental ν^expt_i normal mode vibrational frequencies of methane (in cm⁻¹). Activity in infrared (IR), Raman (R) and hyper-Raman (hR) spectroscopies given for T_d symmetry. Symmetry character Γ and degeneracy n_deg are also given for each mode

Mode	Γ	n deg	ω i	ν ^theo i	ν ^expt i	Act. in
a Experimental data from ref. 66.
1	A 1	1	3026	2903	2917	R, hR
2	E	2	1556	1521	1534	R
3	F 2	3	3126	2990	3019	IR, R, hR
4	F 2	3	1339	1307	1306	IR, R, hR

---

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.

Fig. 1 Vibrational spectra of methane. (a) Theoretical (top) infrared and (bottom) non-polarized Raman spectrum. (b) Theoretical non-polarized hyper-Raman spectrum. (c) Experimental non-polarized hyper-Raman spectrum from ref. 65.

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¹ that is expected at 1339 cm⁻¹ 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¹, 3¹ and 4¹. The 3¹ 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.³⁹ However, including anharmonic corrections increases the difference in intensity between the strong 3¹ peak and the two other weak ones. In the experimental hyper-Raman spectrum, one feature appears larger than the experimental noise around 3000 cm⁻¹. This feature corresponds to the intense 3¹ peak in the theoretical spectra. A shoulder on the lower-frequency side of the feature (around 2900 cm⁻¹) appears in the experimental spectrum, and can be linked to the 1¹ peak predicted in the theoretical spectra. In the experimental reference article,⁶⁵ a peak is also mentioned at 1306 cm⁻¹, corresponding to 4¹. 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¹ and 4¹ 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¹ peak.

4.2 Ethane

The calculated harmonic and anharmonic vibrational frequencies for ethane are compiled in Table 2 and compared to experimental values.⁶⁶ 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 pseudo-degeneracy between modes 1¹ and 5¹ (3024 cm⁻¹), both of these modes corresponding to symmetric stretching motion of the CH₃ moieties (one fully symmetric, the other antisymmetric with respect to the inversion). The anharmonic treatment removes this accidental degeneracy (1¹ at 2868 cm⁻¹ and 5¹ at 2877 cm⁻¹). While experimentally 1¹ is higher in frequency than 5¹, 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¹, 3¹ and 12¹ with absolute deviations from the experimental value going from 70 to 86 cm⁻¹, from 1 to 24 cm⁻¹ and from 3 to 4 cm⁻¹, respectively. Apart from these three modes, the harmonic treatment shows an overestimation of the vibrational frequencies compared to experiment of more than 90 cm⁻¹ for the modes above 2000 cm⁻¹, between 30 and 35 cm⁻¹ for the modes in the 1000 to 2000 cm⁻¹ frequency range, and less than 15 cm⁻¹ for the lower-frequency modes. Adding anharmonic corrections reduces the absolute deviation from experiment to less than 40 cm⁻¹ for the modes above 2000 cm⁻¹ and less than 10 cm⁻¹ for the modes vibrating at lower frequencies.

Table 2 Theoretical harmonic ω_i, anharmonic ν^theo_i and experimental ν^expt_i normal mode vibrational frequencies of ethane (in cm⁻¹). Activity in infrared (IR), Raman (R) and hyper-Raman (hR) spectroscopies given for D_3d symmetry. Symmetry character Γ and degeneracy n_deg are also given for each mode

Mode	Γ	n deg	ω i	ν ^theo i	ν ^expt i	Act. in
a Experimental data from ref. 66.
1	A 1g	1	3024	2868	2954	R
2	A 1g	1	1420	1393	1388	R
3	A 1g	1	994	971	995	R
4	A 1u	1	302	279	289	hR
5	A 2u	1	3024	2877	2896	IR, hR
6	A 2u	1	1411	1382	1379	IR, hR
7	E g	2	3066	2932	2969	R
8	E g	2	1501	1457	1468	R
9	E g	2	1221	1193	1190	R
10	E u	2	3091	2955	2985	IR, hR
11	E u	2	1504	1463	1469	IR, hR
12	E u	2	825	826	822	IR, hR

---

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.

Fig. 2 Vibrational spectra of ethane. (a) Theoretical (top) infrared and (bottom) non-polarized Raman spectrum. (b) Theoretical non-polarized hyper-Raman spectrum. (c) Experimental non-polarized hyper-Raman spectrum from ref. 65.

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¹ at 1411 cm⁻¹ for IR and 2¹ at 1420 cm⁻¹ and 9¹ at 1221 cm⁻¹ 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¹ peak at 302 cm⁻¹ 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¹ 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⁻¹, corresponding to the 8¹11¹ combination band. Two very weak peaks can also be seen at 2745 and 3755 cm⁻¹ corresponding, respectively, to the 2¹6¹ and 7¹12¹ combination modes. In the experimental spectrum, four features can be distinguished: a strong band around 3000 cm⁻¹ and three weak ones that are hardly discernible from the experimental noise around 1400, 800 and 300 cm⁻¹. The main feature corresponds to the group of peaks predicted by the theoretical approach in this region (5¹, 10¹ and 8¹11¹). The three other peaks correspond respectively to the two-peak group 6¹ and 11¹, to the 12¹ peak and to the 4¹ peak. The global intensity of the group of peaks around 3000 cm⁻¹ 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 lower-intensity features are similar to each other.

4.3 Ethylene

The calculated harmonic and anharmonic vibrational frequencies for ethylene are compiled in Table 3 and compared to experimental values.⁶⁶ 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¹ peak at 1245 cm⁻¹, is also predicted but with a relative intensity so low that it does not appear in the spectrum. In hyper-Raman, the 12¹ peak at 1478 cm⁻¹ also has a very low intensity and is almost not visible in the theoretical spectrum. The 4¹ peak at 1060 cm⁻¹ is present in the hyper-Raman spectrum, but neither in the IR nor in the Raman spectra.

Table 3 Theoretical harmonic ω_i, anharmonic ν^theo_i and experimental ν^expt_i normal mode vibrational frequencies of ethylene (in cm⁻¹). Activity in infrared (IR), Raman (R) and hyper-Raman (hR) spectroscopies given for D_2h symmetry. Symmetry character Γ and degeneracy n_deg are also given for each mode

Mode	Sym.	Deg.	ω i	ν ^theo i	ν ^expt i	Act. in
a Experimental data from ref. 66.
1	A g	1	3137	3004	3026	R
2	A g	1	1688	1645	1623	R
3	A g	1	1380	1355	1342	R
4	A u	1	1060	1051	1023	hR
5	B 1g	1	3192	3052	3103	R
6	B 1g	1	1245	1226	1236	R
7	B 1u	1	984	970	949	IR, hR
8	B 2g	1	976	963	943	R
9	B 2u	1	3220	3077	3106	IR, hR
10	B 2u	1	834	835	826	IR, hR
11	B 3u	1	3123	2973	2989	IR, hR
12	B 3u	1	1478	1444	1444	IR, hR

---

The harmonic treatment overestimates vibrational frequencies by more than 30 cm⁻¹ for all modes but 6¹ and 10¹, which are overestimated by less than 10 cm⁻¹. For the modes with the higher frequencies (above 2900 cm⁻¹), the overestimation is larger than 85 cm⁻¹. All anharmonic frequencies are in agreement with the experimental ones, with absolute deviations lower than 30 cm⁻¹, with the exception of mode 5¹ for which the deviation is 51 cm⁻¹. Only the two normal modes with harmonic deviation less than 10 cm⁻¹ (6¹ and 10¹) 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⁻¹, still remaining lower than 10 cm⁻¹.

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.

Fig. 3 Vibrational spectra of ethylene. (a) Theoretical (top) infrared and (bottom) non-polarized Raman spectrum. (b) Theoretical non-polarized hyper-Raman spectrum. (c) Experimental non-polarized hyper-Raman spectrum from ref. 65.

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⁻¹ arising from the 4¹7¹ combination mode. In the experimental spectrum, two main features can be observed around 1000 cm⁻¹ and around 3000 cm⁻¹. These two features can be linked respectively to the 4¹, 7¹ and 10¹ three-peak group and the 9¹ and 11¹ two-peak group. The intensity of the 1000 cm⁻¹ group is higher than the one of the 3000 cm⁻¹ group in a similar way that the 1000 cm⁻¹ feature is more intense than the 3000 cm⁻¹ 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.

Fig. 4 Hyper-Raman spectra of ethylene in HV polarization. (a) Theoretical spectrum. (b) Experimental spectrum from ref. 65.

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¹6¹ combination mode at 2608 cm⁻¹ and the 5¹7¹ combination mode at 3962 cm⁻¹. Also, the relative intensities of the two groups of peaks around 1000 cm⁻¹ and 3000 cm⁻¹ in the theoretical spectra are in agreement with the experimental ones, but the dominant contributing mode to the feature around 1000 cm⁻¹ changes from 7¹ to 4¹ 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 highest-order 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.

5 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 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¹11¹ combination peak at 2914 cm⁻¹ 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.⁶⁹

Acknowledgements

This work has received support from the Research Council of Norway through a Centre of Excellence Grant (Grant No. 179568/V30) and a research grant (Grant No. 191251), the European Research Council through an ERC Starting Grant (Grant No. 279619), and from the Norwegian Supercomputer Program through a grant of computer time (Grant No. NN4654K).

References

1. J. C. Decius and J. E. Rauch, at the Ohio State Symposium on molecular Structure and Spectroscopy, paper 48, 1959.
2. R. W. Terhune, P. P. Maker and C. M. Savage, Phys. Rev. Lett., 1965, 14, 681.
3. A. M. Kelley, Annu. Rev. Phys. Chem., 2010, 61, 41.
4. W. P. Acker, D. H. Leach and R. K. Chang, Chem. Phys. Lett., 1989, 155, 491.
5. A. Mohammed, H. Ågren, M. Ringholm, A. J. Thorvaldsen and K. Ruud, Mol. Phys., 2012, 110, 2315.
6. S. Bratoz, Colloq. Int. C. N. R. S., 1958, 82, 287.
7. P. Pulay, Mol. Phys., 1969, 17, 197.
8. J. A. Pople, R. Krishnan, H. B. Schlegel and J. S. Binkley, Int. J. Quantum Chem., Quantum Chem. Symp., 1979, S13, 225.
9. B. Brooks, W. D. Laidig, P. Saxe, J. Goddard, Y. Yamaguchi and H. F. Schaefer, J. Chem. Phys., 1980, 72, 4652.
10. R. Krishnan, H. B. Schlegel and J. A. Pople, J. Chem. Phys., 1980, 72, 4654–4655.
11. M. R. Hoffmann, D. J. Fox, J. F. Gaw, Y. Osamura, Y. Yamaguchi, R. S. Grev, G. B. Fitzgerald, H. F. Schaefer, P. J. Knowles and N. C. Handy, J. Chem. Phys., 1984, 80, 2660–2668.
12. P. Jørgensen and J. Simons, J. Chem. Phys., 1983, 79, 334–357.
13. T. U. Helgaker, J. Almlöf, H. J. Aa. Jensen and P. Jørgensen, J. Chem. Phys., 1986, 84, 6266–6279.
14. N. C. Handy, R. D. Amos, J. F. Gaw, J. E. Rice and E. D. Simandiras, Chem. Phys. Lett., 1985, 120, 151–158.
15. R. J. Harrison, G. B. Fitzgerald, W. D. Laidig and R. J. Bartlett, Chem. Phys. Lett., 1986, 124, 291–294.
16. H. Koch, H. J. Aa. Jensen, P. Jørgensen, T. Helgaker, G. E. Scuseria and H. F. Schaefer, J. Chem. Phys., 1990, 92, 4924–4940.
17. J. Gauss and J. F. Stanton, Chem. Phys. Lett., 1997, 276, 70–77.
18. P. G. Szalay, J. Gauss and J. F. Stanton, Theor. Chem. Acc., 1998, 100, 5–11.
19. A. Komornicki and G. B. Fitzgerald, J. Chem. Phys., 1993, 98, 1398–1421.
20. B. G. Johnson and M. J. Frisch, Chem. Phys. Lett., 1993, 216, 133–140.
21. P. Deglmann, F. Furche and R. Ahlrichs, Chem. Phys. Lett., 2002, 362, 511–518.
22. G. Placzek, Handbuch der Radiologie, 1934, 6(2), 205.
23. E. B. Wilson Jr., J. C. Decius and P. C. Cross, Molecular Vibrations, General Publishing Company, Ltd., Don Mills, Toronto, Canada, 1955.
24. O. Quinet and B. Champagne, J. Chem. Phys., 2003, 117, 2481.
25. O. Quinet and B. Champagne, Theor. Chem. Acc., 2004, 111, 390.
26. O. Quinet, B. Champagne and V. Rodriguez, J. Chem. Phys., 2004, 121, 4705.
27. O. Quinet, B. Champagne and S. J. A. van Gisbergen, Int. J. Quantum Chem., 2006, 106, 599.
28. O. Quinet, B. Champagne and V. Rodriguez, J. Chem. Phys., 2006, 124, 244312.
29. M. Ringholm, R. Bast, L. Oggioni, U. Ekström and K. Ruud, J. Chem. Phys., 2014, 141, 134107.
30. H. H. Nielsen, Rev. Mod. Phys., 1951, 23, 90.
31. D. Papoušek and M. R. Aliev, Molecular Vibrational-Rotational Spectra, Elsevier, Amsterdam, NL, 1981.
32. J. M. L. Martin, T. J. Lee, P. R. Taylor and J.-P. François, J. Chem. Phys., 1995, 103, 2589.
33. J. Neugebauer and B. A. Hess, J. Chem. Phys., 2003, 118, 7215.
34. J. Bloino and V. Barone, J. Chem. Phys., 2012, 136, 124108.
35. J. Breidung, W. Thiel, J. Gauss and J. F. Stanton, J. Chem. Phys., 1999, 110, 3687.
36. J. F. Stanton and J. Gauss, Int. Rev. Phys. Chem., 2000, 19, 61–95.
37. P. Carbonniere and V. Barone, Chem. Phys. Lett., 2004, 399, 226.
38. V. Barone, J. Chem. Phys., 2005, 122, 014108.
39. Y. Cornaton, M. Ringholm, O. Louant and K. Ruud, Phys. Chem. Chem. Phys., 2016, 18, 4201.
40. T.-C. Jagau, J. Gauss and K. Ruud, J. Chem. Phys., 2013, 139, 154106.
41. P. E. Maslen, D. Jayatilaka, S. M. Colwell, R. D. Amos and N. C. Handy, J. Chem. Phys., 1991, 95, 7409.
42. P. E. Maslen, N. C. Handy, R. D. Amos and D. Jayatilaka, J. Chem. Phys., 1992, 97, 4233.
43. M. Ringholm, D. Jonsson and K. Ruud, J. Comput. Chem., 2014, 35, 622.
44. A. J. Thorvaldsen, K. Ruud, K. Kristensen, P. Jørgensen and S. Coriani, J. Chem. Phys., 2008, 129, 214108.
45. B. Gao, A. J. Thorvaldsen and K. Ruud, Int. J. Quantum Chem., 2011, 111, 858.
46. S. Reine, E. I. Tellgren and T. Helgaker, Phys. Chem. Chem. Phys., 2011, 9, 4771.
47. U. Ekström, L. Visscher, R. Bast, A. J. Thorvaldsen and K. Ruud, J. Chem. Theory Comput., 2010, 6, 1971.
48. M. Ringholm, D. Jonsson, R. Bast, B. Gao, A. J. Thorvaldsen, U. Ekström, T. Helgaker and K. Ruud, J. Chem. Phys., 2014, 140, 034103.
49. S. J. Cyvin, J. E. Rauch and J. C. Decius, J. Chem. Phys., 1965, 43, 4083.
50. V. Barone, M. Biczysko and J. Bloino, Phys. Chem. Chem. Phys., 2014, 16, 1759.
51. K. Kristensen, P. Jørgensen, A. J. Thorvaldsen and T. Helgaker, J. Chem. Phys., 2008, 129, 214103.
52. R. Bast, XCINT library, http://dftlibs.org/xcint/, 2016.
53. K. Aidas, C. Angeli, K. L. Bak, V. Bakken, R. Bast, L. Boman, O. Christiansen, R. Cimiraglia, S. Coriani, P. Dahle, E. K. Dalskov, U. Ekström, T. Enevoldsen, J. J. Eriksen, P. Ettenhuber, B. Fernández, L. Ferrighi, H. Fliegl, L. Frediani, K. Hald, A. Halkier, C. Hättig, H. Heiberg, T. Helgaker, A. C. Hennum, H. Hettema, E. Hjertenæs, S. Høst, I.-M. Høyvik, M. F. Iozzi, B. Jansík, H. J. A. Jensen, D. Jonsson, P. Jørgensen, J. Kauczor, S. Kirpekar, T. Kjærgaard, W. Klopper, S. Knecht, R. Kobayashi, H. Koch, J. Kongsted, A. Krapp, K. Kristensen, A. Ligabue, O. B. Lutnæs, J. I. Melo, K. V. Mikkelsen, R. H. Myhre, C. Neiss, C. B. Nielsen, P. Norman, J. Olsen, J. M. H. Olsen, A. Osted, M. J. Packer, F. Pawlowski, T. B. Pedersen, P. F. Provasi, S. Reine, Z. Rinkevicius, T. A. Ruden, K. Ruud, V. V. Rybkin, P. Sałek, C. C. M. Samson, A. S. de Merás, T. Saue, S. P. A. Sauer, B. Schimmelpfennig, K. Sneskov, A. H. Steindal, K. O. Sylvester-Hvid, P. R. Taylor, A. M. Teale, E. I. Tellgren, D. P. Tew, A. J. Thorvaldsen, L. Thøgersen, O. Vahtras, M. A. Watson, D. J. D. Wilson, M. Ziolkowski and H. Ågren, Wiley Interdiscip. Rev.: Comput. Mol. Sci., 2013, 4, 269.
54. DALTON, a molecular electronic structure program. Release DALTON2013, see http://daltonprogram.org/, 2013.
55. P. Jørgensen, H. J. A. Jensen and J. Olsen, J. Chem. Phys., 1988, 89, 3654.
56. B. Gao and A. J. Thorvaldsen, GEN1INT Version 0.2.1, a library to evaluate the derivatives of one-electron integrals with respect to the geometry perturbations, external electric and magneti fields, and total rotational angular momentum at zero fields with contracted rotational London atomic orbitals, released under the GNU Lesser General Public License, 2012.
57. A. J. Thorvaldsen, CGTO-DIFF-ERI, a library for the evaluation of geometry-differentiated two-electron integrals (dⁿ/dRⁿ_k) under the GNU Lesser General Public License, 2012.
58. U. Ekström, XCFUN library, http://dftlibs.org/xcfun/, 2010.
59. J. H. Meal and S. R. Polo, J. Chem. Phys., 1956, 24, 1126.
60. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098.
61. C. Lee, W. Yang and R. G. Paar, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785.
62. P. J. Stephens, F. J. Delvin, C. F. Chabalowski and M. J. Frisch, J. Phys. Chem., 1994, 98, 11623.
63. T. H. Dunning Jr., J. Chem. Phys., 1989, 90, 1007.
64. I. M. Mills, in Molecular Spectroscopy: Modern Research, ed. K. N. Rao and C. W. Mathews, Academic, New York, NY, USA, 1972, pp. 115–140.
65. J. F. Verdieck, S. H. Peterson, C. M. Savage and P. D. Maker, Chem. Phys. Lett., 1970, 7, 219.
66. T. Shimanouchi, Tables of Molecular Vibrational Frequencies Consolidated, National Bureau of Standards, 1972, vol. 1, pp. 1–160.
67. D. Begue, P. Carbonniere, V. Barone and C. Pouchan, Chem. Phys. Lett., 2005, 415, 25.
68. C. Puzzarini, M. Biczysko and V. Barone, J. Chem. Theory Comput., 2010, 6, 828.
69. N. H. List, M. T. P. Beerepoot, J. M. H. Olsen, B. Gao, K. Ruud, H. J. A. Jensen and K. Kongsted, J. Chem. Phys., 2015, 142, 034119.

---