Cumulant expansion for fast estimate of non-Condon effects in vibronic transition profiles.

When existing, cumulants can provide valuable information about a given distribution and can in principle be used to either fully reconstruct or approximate the parent distribution function. A previously reported cumulant expansion approach for Franck-Condon profiles [Faraday Discuss., 150, 363 (2011)] is extended to describe also the profiles of vibronic transitions that are weakly allowed or forbidden in the Franck-Condon approximation (non-Condon profiles). In the harmonic approximation the cumulants of the vibronic profile can be evaluated analytically and numerically with a coherent state-based generating function that accounts for the Duschinsky effect. As illustration, the one-photon 1 1Ag → 1 1B2u UV absorption profile of benzene in the electric dipole and (linear) Herzberg-Teller approximation is presented herein for zero Kelvin and finite temperatures.

Introduction

Vibrationally resolved electronic spectra (e.g. one-photon absorption and emission spectra) are within the Born-Oppenheimer framework usually interpreted in terms of Franck–Condon (FC) factors (FCFs)[1,2]. Accordingly, one can try to obtain the shape of the spectral profile for a FC-allowed transition from computed FCFs in frequency domain. However, besides the electronic structure calculations, the evaluation of FCFs for large molecular systems is challenging even within the harmonic approximation if one has to take Duschinsky mode mixing (rotation)[3] into account. This is because multi-variate Hermite polynomials have then to be evaluated for each FC integral, rather than only uni-variate Hermite polynomials as is the case for the comparatively simple parallel harmonic oscillator model. The computational task becomes more difficult as the molecular size and temperature increases because the number of FC integrals grows vastly. The evaluation of the multi-dimensional FC integral in the harmonic approximation is, at least for some cases with Duschinsky mode mixing, classified as a #P-hard problem in computational complexity theory[4] and thus it recently became a topic in quantum computation. A quantum optical simulation (quantum computation) has been proposed theoretically for instance for the FC profile calculation and has been performed for the photoelectron spectrum of SO2 by a trapped-ion device[5-7].

To describe FC-forbidden or weakly allowed transitions, one has to go beyond the Condon approximation and employ for instance a Herzberg–Teller (HT) expansion[8] of the electronic transition moment with respect to the normal coordinates. As a result, the calculation of the vibronic spectrum for a non-Condon process is even more difficult than for a FC-allowed transition because one has to evaluate matrix elements of the non-Condon operators which require for each HT integral in general the calculation of combinations of several FC integrals.

The number of FC integrals and matrix elements of non-Condon operators to be evaluated in a sum-over-states approach can, in principle, be significantly reduced with the help of integral configuration selection strategies[9,10]. However, this time-independent (TI) calculation of the spectral profile in the frequency domain is still considerably more expensive than an alternative time-dependent (TD) approach that exploits time-correlation functions (TCFs) (see e.g. ref.[11-15]), but offers the ability to directly assign individual peaks in the spectrum. As we have outlined earlier[16], a unified coherent state-based generating function (CSGF) approach[17] can be used both for rigorous integral prescreening strategies and TCF calculations which combine the strengths of both approaches and complement each other favorably. Even in the less demanding TD approach, however, one usually invests significant computational time for often unnecessary spectral details.

Cumulants (or moments) of a distribution (see e.g. refs[17-27]) can deliver highly useful information. From this one can either attempt to reconstruct the spectral shape or try to estimate the relevant spectral profile, which can be exploited in subsequent TI and TD approaches[16,28]. Cumulants of the vibronic spectrum can be obtained from the CSGF directly without computing the total spectrum in frequency domain. This method was exploited already in ref.[29] for FC-allowed transitions, and we report herein an extension of this method to incorporate non-Condon transitions. To illustrate the performance of the approach, we present the profile of the 1A1g → 1B2u transition of benzene, which frequently served as a prototypical example for multiple authors (see e.g. refs[30-32] and references therein). The transition is in the electric dipole approximation Franck–Condon forbidden and it is studied herein at various temperatures within the linear HT and harmonic approximation. Cumulant expansion is compared herein to the TCF approach.

Results

Theory

The spectral profile (ρ(ħ ω; T)) can be expressed via the Fourier transform (FT) of the TCF (χ(t; T)) that depends on the time t and temperature T, namelywhere ħω is the transition energy and ħω 0 the 0′ − 0 transition energy. The corresponding occupancy representation for the TCF and the spectral profile can be obtained from Fermi’s Golden Rule, respectively, asandwhere we have assumed an electric dipole transition with the electronic transition dipole moment , which is a function of normal coordinates of the initial electronic state, and the harmonic approximation. The N-dimensional harmonic oscillator eigenstates of the initial and final electronic state are denoted by and with the corresponding harmonic energy vectors and , respectively. and are the N-dimensional harmonic oscillator Hamiltonians belonging to the initial and final electronic states, respectively. k B is the Boltzmann constant. is the vibronic transition energy with respect to the 0′ − 0 transition energy. The spatial representation of the TCF in closed form can be found by evaluating the following quantum mechanical traces

The traces can be evaluated with any complete basis. In our work N-dimensional coherent states were used (see e.g. refs[9,16]) with the Duschinsky relation between initial and final state normal coordinates ( where and are the Duschinsky rotation matrix and displacement vector, respectively and are the normal coordinates of the final state) and the linear HT expansion of the electronic transition dipole moment ( where is the first derivative of with respect to .).

The TCF is related to a probability density function (PDF). If all cumulants or moments of a PDF are defined and available, the PDF can be reconstructed as follows[29] where is the k-th order cumulant at temperature T. The cumulants of the spectral density function are normalised to the total integrated profile [16]. Moments (cumulants and moments are inter-convertible[33]) can be obtained by partial derivatives of χ with respect to time,

The cumulants can be obtained from the moments via the following transformation[33],

Thus, cumulants can be evaluated analytically or numerically by evaluating partial derivatives of χ in Eq. (4) with respect to the time variable at t = 0. Analytic evaluation of the cumulants to arbitrary order within the linear HT approximation can be performed along the lines of the development in refs[16,28,29] for the cumulants of FC profiles to arbitrary order. For numerical evaluation of low-order cumulants one needs to compute χ at the first few time steps. To obtain the corresponding moments numerically, Re(χ(t, T)) and Im(χ(t, T)) as computed at these time steps are used to determine low-order even and odd moments, respectively, because in Eq. (2) (see e.g. ref.[23]).

The closed form of χ(t, T) within the linear HT approximation can be found in refs[14-16,28]. In the present work, flexibility is used in the GF to obtain detailed information concerning individual contributions of different modes. This is achieved by assigning different time and temperature variables to each vibrational degree of freedom. The corresponding GF in an occupancy representation reads as followswhere the general operators and , which can be products of momentum and position operators, are given instead of ; and different temperatures can be given to the initial and final vibrational degrees of freedom viawhere β  = 1/(k B T ) (with Boltzmann constant k B and temperature T ). The parameter matrices are defined as followswith the time variables being assigned to the matrices and for initial and final vibrational modes, respectively, andwith .

is the corresponding normalizing factor related to the partition function of the Boltzmann distribution of harmonic oscillators, i.e.

The Duschinsky relation is considered with the Doktorov matrices and vectors[34] as well as

The temperature dependent parameters are defined as follow

Finally, the Franck–Condon Herzberg–Teller (FCHT) TCF reads aswith the FC generating function,the mixed FC/HT generating function,and the HT generating function,

The electronic 1 1Ag → 1 1B2u transition of benzene is FC-forbidden in the electric dipole approximation such that only the HT terms contribute to the spectral function. The corresponding TCF for FCHT weighted density of states (FCHTW) is given as follows, here with same time (t) for all vibrational degrees of freedom and same temperature (T  = T and T ′ = ∞) for the initial and final modes, respectively,

Accordingly, the spectrum is obtained by the Fourier transformation,which can show the detailed vibronic structure.

The FCHTW profile is now approximated with a finite number of cumulants via the Edgeworth expansion with the order n [35]. Whereas for n = 2 a Gaussian distribution function is used, the Edgeworth expansion for order n ≥ 3 is employed as[29] where is a set of non-negative integer vectors, which are constrained to and . is a uni-variate Hermite polynomial of order s + 2r and is defined as follows,

The Edgeworth expansion with a finite number of cumulants and in the infinite series are related asand we use the relation, , for the wavenumber domain profile.

Computational detail

The vibronic profiles for benzene’s 1 1Ag → 1 1B2u transition at zero Kelvin and finite temperatures are calculated with the two methods, namely TCF and time-independent cumulant expansion (CE). We use herein the term time-independent CE, which was employed in ref.[29], to distinguish this CE from the conventional (time-dependent) CE (see e.g. refs[20-24]) which involves time integration for the cumulant calculation. To compute the vibronic spectra via the TCF method, the FFTW library[36] is used for fast Fourier Transform (FFT). The approximate curves are generated for the CE with Edgeworth expansion[29,35] using the computed low-order cumulants. Some of the problems related to this type of expansion for the description of FC profiles are discussed in ref.[29]. The moments (Eq. 6) are calculated both analytically and numerically for comparison, the latter by approximating partial derivatives of χ in Eq. 4 with respect to time (see results section) via a central finite difference scheme with a truncation error being of order (Δt)2. When generating the data points in time, we exploit the time-reversal symmetry condition, i.e. χ(−t; T) = χ(t; T)*. Required input data from electronic structure calculations for benzene, i.e. molecular equilibrium structures and corresponding harmonic force fields for each electronic state (1Ag and 1B2u) as well as first derivatives of the electronic transition dipole moments are taken from ref.[30] (CASSCF/DZV). These data have been compared to results obtained via analytical derivative techniques for electronic transition dipole moments within a time-dependent density functional theory framework in ref.[32]. The vibronic structure methods employed in the present work are implemented in a development version of our vibronic structure program package hotFCHT[9,30].

Numerical simulation

The computed vibronic spectra are shown in Fig. 1. The left hand side in Fig. 1 shows vibronic profiles from TCF-FFT which are convoluted by a Lorentzian line shape function with full width at half maximum (FWHM) of 50 cm−1 at temperatures elevating from 0 K to 1000 K. This FC-forbidden vibronic transition is mediated by the non-totally symmetric vibrational modes in the irreducible representation e2g of the D6h molecular symmetry group. The main feature of the vibronic spectrum is from progressions in the totally symmetric C-C stretching mode (963 cm−1) building on the so-called false origin from a single excitation of a non-totally symmetric (e2g) in-plane bending mode (575 cm−1) as indicated in the spectrum at zero Kelvin. The calculated spectrum at 300 K is compared with the experimental data of Fischer[37]. The two spectra agree fairly well in the low energy region but the computed peaks at higher energies are slightly shifted to larger wavenumbers due to the harmonic approximation. As temperature increases the vibrational structure becomes very congested and washed out. At 1000 K (only employed for testing the method) one can not see a resolved vibrational structure any longer, only the corresponding envelope.

Figure 1

Left part of the figure: The dashed lines are drawn for the TCF-FFT approach with a Lorentzian line shape function with FWHM of 50 cm−1. A time increment Δt of 0.51 fs and a grid with 216 grid points are used for the corresponding FFT calculations. The experimental UV absorption spectrum as reported by Fischer in ref.[37] is additionally shown in red, which has been adapted from ref.[37] and shifted to match approximately the position of the major peak in the region below the 0–0 transition wavenumber and rescaled to have similar peak height as the one computed for the transition. Right part of the figure: The dashed lines are drawn for the TCF-FFT approach with a Gaussian line shape function of with FWHM of 500 cm−1. A time increment of 0.10 fs and a grid with 215 grid points are used for the corresponding FFT calculations. Solid lines are drawn for the curve obtained by Edgeworth expansion using up to 4-th order cumulants and, for 1000 K by Edgeworth expansion using up to 8-th order cumulants.

On the right hand side of Fig. 1 the two methods (TCF-FFT and CE-Edgeworth) are compared for increasing temperatures. The spectra are convoluted in the TCF-FFT curves (dashed lines) with a Gaussian line shape function of 500 cm−1 for FWHM. The second moment (4.51 × 104 cm−2 (hc 0)2) of the Gaussian line shape function is added to the second moments of the vibronic spectrum to take the line shape function into account (see ref.[29] for the rationalisation and details). The relatively broad line shape function is used for the TCF-FFT curves to have vibrationally relatively structureless spectra for comparison. At 0, 300 and 500 K, the TCF-FFT curves still show vibrational structure and the CE-Edgeworth curves (solid lines) look like nonlinear regression curves of the corresponding TCF-FFT spectra. When the vibrational structures are also essentially smoothed out in the TCF-FFT curves at 1000 K, the two approaches agree with each other extremely well. Up to the 4-th order cumulants are used for 0, 300, 500 K and up to 8-th order cumulants are computed for 1000 K.

In Table 1 the moments computed numerically (via numerical derivatives) and analytically are compared. At low orders and all temperatures the two methods agree well and for higher orders still the agreement is satisfactory. One of the advantages of the numerical method is that one only needs to compute the TCF for the first few time steps and it can be improved by controlling the time increment and the number of data points. The analytical method usually meets a combinatorial problem in high order cumulant calculations due to the analytic derivatives of the inverse matrix[29]. The second advantage of the numerical method is that it is easy to include linear and nonlinear non-Condon effects. The third advantage is that one can incorporate general line shape functions which would not have well defined cumulants (see the discussion on page 415 of ref.[29]). Lastly, the computational cost of the numerical cumulant expansion method is that the number of data points to be evaluated is almost negligible comparing to the TCF-FFT approach, which is about three orders of magnitude more expensive.

Table 1

Analytically and numerically computed moments.

n	〈Eε_′,ε_n〉/(cm^−1hc_0)^n
	T = 0 K		T = 300 K		T = 500 K		T = 1000 K
	Analytical	Numerical	Analytical	Numerical	Analytical	Numerical	Analytical	Numerical
1	2.61 × 103	2.61 × 103	2.47 × 103	2.47 × 103	2.12 × 103	2.12 × 103	1.12 × 103	1.12 × 103
2	9.24 × 106	9.21 × 106	8.64 × 106	8.62 × 106	7.38 × 106	7.36 × 106	5.64 × 106	5.63 × 106
3	4.07 × 1010	4.05 × 1010	3.77 × 1010	3.75 × 1010	3.17 × 1010	3.15 × 1010	2.03 × 1010	2.01 × 1010
4	2.14 × 1014	2.12 × 1014	1.97 × 1014	1.95 × 1014	1.66 × 1014	1.64 × 1014	1.29 × 1014	1.27 × 1014
5	—	—	—	—	—	—	8.00 × 1017	7.82 × 1017
6	—	—	—	—	—	—	6.55 × 1021	6.33 × 1021
7	—	—	—	—	—	—	5.84 × 1025	5.54 × 1025
8	—	—	—	—	—	—	6.16 × 1029	5.71 × 1029

4.51 × 104 cm−2 (hc 0)2 is added to the second moments to take the Gaussian line shape function (FWHM = 500 cm−1) into account; see ref.[29] for details. A time increment of 0.10 fs is used for computing the numerical derivatives.

Mean excitation wavenumbers ( with being a number operator of i-th mode in the final electronic state) of individual modes are computed analytically for the HT active e2g symmetric vibrational modes of the final state, and are given in Table 2. The corresponding first derivative in Eq. (6) can be performed numerically or analytically by assigning z = diag(1, …, 1) and to Eq. 17. The mean excitation energy can serve as a parameter for the individual vibrational degrees of freedom as an effective reorganisation energy or a Huang–Rhys factor (when normalised by its harmonic energy), which can be characterised as a function of structural deformation, frequency change, Duschinsky mode coupling and temperature both in the Condon and non-Condon approximation. One might naively expect a larger mean excitation energy as the temperature increases but the mean values of high frequency modes (1665 and 3389 cm−1) in some intermediate temperature ranges are smaller than at zero Kelvin. This can be rationalised as follows: Because the Duschinsky mode mixing between low and high frequency modes is small in the present case, the high frequency modes can not obtain thermal energy from the low frequency modes efficiently. Thus the high frequency modes are almost thermally inactive, whereas the total integrated profile (ρ tot) increases as temperature increases. In the mean energy calculation of the high frequency modes at finite temperatures the denominator (total intensity) increases because low frequency modes accept thermal energy while the numerator (excitation of high frequency modes) stays constant. Therefore the mean excitation energies of high frequency modes are reduced at finite temperatures. If Duschinsky rotation couples the low and high frequency modes significantly, however, thermal energy can be transfered to the high frequency modes via the low frequency modes in the initial state, accordingly the mean excitation energies of high frequency modes can increase as temperature increases.

Table 2

Mean excitation wavenumbers of the components of individual vibrational e2g symmetric modes of benzene as computed for different temperatures.

ε_i′〈vˆ_i′〉/(hc_0cm^−1)
Mode ν i	ω˜_e′/cm−1	T = 0 K	T = 300 K	T = 500 K	T = 1000 K
ν 6	575	2.62 × 102	2.90 × 102	3.79 × 102	7.25 × 102
ν 6	575	2.62 × 102	2.90 × 102	3.79 × 102	7.25 × 102
ν 9	1237	1.60 × 101	1.78 × 101	4.64 × 101	2.51 × 102
ν 9	1237	1.60 × 101	1.78 × 101	4.64 × 101	2.51 × 102
ν 8	1665	2.54 × 101	2.38 × 101	3.09 × 101	1.65 × 102
ν 8	1665	2.54 × 101	2.38 × 101	3.09 × 101	1.65 × 102
ν 7	3389	8.14 × 101	7.49 × 101	6.08 × 101	6.37 × 101
ν 7	3389	8.14 × 101	7.49 × 101	6.08 × 101	6.37 × 101

The numbering used for the modes ν 6, ν 7, ν 8 and ν 9, correspond to that used by Wilson for benzene and translates to ν 18, ν 15, ν 16 and ν 17 in Herzberg’s nomenclature, respectively. The corresponding harmonic vibrational wavenumbers as computed in ref.[30] for the electronically excited state and as used in the present calculations are also given.

In closing the section, the moments or the cumulants of the vibronic excitation energy can provide intuitively useful information concerning the vibronic transition profile with almost no computation cost comparing to the TCF-FFT method. Furthermore, the mean excitation energy of individual mode opens a new interpretation for the vibronic transition with a single quantity incorporating the mode mixing and the non-Condon effects as well as the temperature, geometrical change and distortion effects.

Conclusion and outlook

We have discussed a cumulant expansion method for describing non-Condon transitions and applied it to the prototypical one-photon electric dipole 1 1Ag → 1 1B2u transition of benzene, which is FC forbidden but HT allowed in the linear HT approximation. The method is particularly powerful when one does not require all the details of the vibronic structures, but rather only quantities such as peak maximum, mean and variance of the spectral shape. This method is computationally much cheaper than the sum-over-states and time-correlation function approach. Moreover, the information (e.g. on the spectroscopically relevant energy window) from the cumulant expansion method can be used in the calculation within the other two methods. For example, the relevant vibronic transition energy domain obtained from the cumulant expansion method can be used for screening the excitation configurations to avoid the unnecessary integral calculations in the time-independent approach. In the time-dependent approach, the same information can be exploited to define via the reciprocal energy window the time increment for the numerical propagation of the wavepacket. The number of discrete time steps with a given time increment determines the desired spectral resolution and the total wavepacket propagation time. Herein we compared a numerical approach for the calculation of cumulants with the results from an analytical scheme. The results obtained numerically are still fairly good. With this method, one can incorporate easily nonlinear non-Condon terms and various line shape functions. In the time-correlation function calculation the real part and imaginary part at each time step provide automatically the even and odd moments, respectively. In the first few time steps we already have the first few moments available and the probability distribution function (information) becomes complete as time progresses. The benzene example selected herein serves to illustrate the principle of the method for future routine applications for molecular systems with hundreds of atoms. As an outlook for the approach described herein, we would like to indicate that it can be generalized for the even more challenging anharmonic vibronic transition problem, because the method has a clear connection to the time correlation function approach. When the time correlation function can be evaluated only for the first few time steps, the cumulants can be determined numerically. The overall computational cost is only the evaluation of a few time correlation function.