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.

Cumulants (or moments) of a distribution (see e.g. refs [17][18][19][20][21][22][23][24][25][26][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 1 A 1g → 1 B 2u transition of benzene, which frequently served as a prototypical example for multiple authors (see e.g. refs [30][31][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, namely where ħω 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, as where we have assumed an electric dipole transition with the electronic transition dipole moment μ Q ( ( )), 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 , respectively. Ĥ and ′ H are the N-dimensional harmonic oscillator Hamiltonians belonging to the initial and final electronic states, respectively. k B is the Boltzmann constant. ε ε ′ E , 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 ( ′ = + Q Q d S where S and d are the Duschinsky rotation matrix and displacement vector, respectively and ′ Q are the normal coordinates of the final state) and the linear HT expansion of the electronic transition dipole moment is the first derivative of µ with respect to Q k .). 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 is the k-th order cumulant at temperature T. The cumulants of the spectral density function are normalised to the total integrated profile ρ 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][15][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 follows where the general operators f and ĝ , which can be products of momentum and position operators, are given instead of μQ ( ); and different temperatures can be given to the initial and final vibrational degrees of freedom via The parameter matrices are defined as follows with the time variables being assigned to the matrices  for initial and final vibrational modes, respectively, and .  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 The temperature dependent parameters are defined as follow T T

Finally, the Franck-Condon Herzberg-Teller (FCHT) TCF reads as
the mixed FC/HT generating function, The electronic 1 1 A g → 1 1 B 2u transition of benzene is FC-forbidden in the electric dipole approximation µ = ( (0) 0) 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 k = T and T k ′ = ∞) 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 is a uni-variate Hermite polynomial of order s + 2r and  + m 2 is defined as follows, Computational detail. The vibronic profiles for benzene's 1 1 A g → 1 1 B 2u 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][21][22][23][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 ( 1 A g and 1 B 2u ) 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 e 2g of the D 6h 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 (e 2g ) 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.
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 × 10 4 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.
Mean excitation wavenumbers (ε ′〈 ′〉 v hc /( ) i i 0 with ′ v i being a number operator of i-th mode in the final electronic state) of individual modes are computed analytically for the HT active e 2g 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 ′ = … … ε ′ z diag(1, , 1, e ,1, 1) 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.
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 1 A g → 1 1 B 2u 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  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 6 0 1 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 2 15 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. 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.   Table 1. Analytically and numerically computed moments. 4.51 × 10 4 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.  8.14 × 10 1 7.49 × 10 1 6.08 × 10 1 6.37 × 10 1 ν 7 3389 8.14 × 10 1 7.49 × 10 1 6.08 × 10 1 6.37 × 10 1 Table 2. Mean excitation wavenumbers of the components of individual vibrational e 2g symmetric modes of benzene as computed for different temperatures. 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 ω ′  e as computed in ref. 30 for the electronically excited state and as used in the present calculations are also given.