Franck-Condon spectra of unbound and imaginary-frequency vibrations via correlation functions: a branch-cut free, numerically stable derivation

For a zero-temperature spectrum, we are interested in calculating the vacuum auto-correlation function,

where $|0\rangle$ is the ground vibrational state of the initial state Hamiltonian, $H_{0}$. (The finite-temperature case will be considered in Section II.2.) The FC spectral density is related to the auto-correlation function by a Fourier transform,

where $\omega$ equals the transition frequency from the $|0\rangle$ initial state, which has zero-point energy $E_{0}=\sum_{k}\omega_{k}/2$.

Our approach to calculating $C_{0}$ is based on disentangling the exponential propagator. First, we recognize that a harmonic Hamiltonian can always be cast as a general quadratic polynomial in the raising and lowering operators, $a_{i}^{\dagger}$ and $a_{i}$, i.e. as

where repeated indices are implied summations and the $n\times n$ matrices $\mathbf{B}$ and $\mathbf{C}$ are symmetric. (For Hermitian $H$, furthermore, $\mathbf{B}=\mathbf{C}$ and $\mathbf{A}$ is symmetric.) The space of quadratic operators is closed under commutation,

thus forming a Lie algebra. This guarantees that exponential quadratic operators can, in principle, be disentangled via the Baker–Campbell–Hausdorff (BCH) formula as

The order of the terms in the argument of the exponent in the first line of this equation is arbitrary, but the order of the individual factors in the second line, Eq. 5, is non-trivial because they do not commute. Assuming the parameters of this particular disentangled form are known, then $C_{0}$ is simple to evaluate,

where we have made repeated use of the fact that $a_{i}|0\rangle=\langle 0|a_{i}^{\dagger}=0$. The problem is now to determine the (time-dependent) disentangled coefficient $h^{\prime}$ and the determinant of $\exp[\mathbf{A}^{\prime}/2]$ given the original entangled exponential operator.

To proceed, we use a faithful matrix representation of the quadratic Lie algebra [25]. Consider the $(2n+2)\times(2n+2)$ matrix, $M$, which is constructed from the quadratic operator coefficients as

It is straightforward to show that this matrix map of the quadratic operators satisfies the same commutation rules. That is, if we let $M$, $M^{\prime}$, and $M^{\prime\prime}$ be the matrices mapped to by the operators $Q$, $Q^{\prime}$, and $Q^{\prime\prime}$, then $[M,M^{\prime}]=M^{\prime\prime}$ if and only if $[Q,Q^{\prime}]=Q^{\prime\prime}$. The faithful matrix representation is useful because BCH formulas can be derived with them via direct matrix exponentiation and multiplication rather than manipulating the exponential operators themselves [26]. It turns out that the matrices of the form $M$ have a relatively straightforward matrix exponential, so that the disentangling rule, Eq. 5, can be expressed as a $(2n+2)\times(2n+2)$ matrix equation. The full derivation is presented in Appendix A. The two main results needed to evaluate the correlation function are

and

where

and

Let us first focus on the matrix $\mathbf{\Gamma}$, which contributes to $h^{\prime}$. The elements in the diagonal positive-frequency matrices $e^{+}$, $\eta^{+}$, and $\zeta^{+}$ that correspond to modes with imaginary frequencies ($\sigma_{k}=-i$) are exponentially divergent as $t\rightarrow\infty$. These divergences, however, ultimately cancel in $\mathbf{\Gamma}$ itself, which we can see by assuming the positive-frequency terms dominate the negative-frequency terms at large time,

The cancellation of exponentially large intermediate quantities implies that the naive evaluation of $\mathbf{\Gamma}$ using finite-precision floating point arithmetic directly as written in Eq. 14 may be unstable and inaccurate. Indeed, numerical tests quickly confirm this issue, suggesting that it is necessary to explicitly remove the divergent intermediates.

To do so, we first inspect various factorizations of the following matrix inverse, which occurs frequently,

where ^′ denotes the transpose-inverse. In each case, the factorization yields a numerically stable matrix inverse and provides another exponentially damped factor on its right or left. After expanding all terms in Eq. 14, one can choose the right- and left-sided forms for each as necessary to cancel the various $\eta^{+}$ and $\zeta^{+}$ factors. After some straightforward, if tedious, algebra, the result is

The non-commuting nature of matrix multiplication makes this expression appear rather cumbersome. Nonetheless, all exponentially divergent intermediates have been removed, and it can be numerically evaluated accurately.

We now proceed to the second factor of $C_{0}(t)$, the determinant of $\exp[\mathbf{A}^{\prime}/2]$. The disentangling formula yields only the matrix exponential of $\mathbf{A}^{\prime}$ itself (Eq. 12). The determinants of these two matrices are, of course, related by a square root, but the relation $\det[e^{\mathbf{A^{\prime}}/2}]=\left(\det[e^{\mathbf{A^{\prime}}}]\right)^{1/2}$ requires careful consideration of the branch-cut choice implicit in the square root in the right-hand side. If $(\cdots)^{1/2}$ is taken to be the principal square root, which we denote with $\sqrt{\cdots}$, then each time $\det[e^{\mathbf{A^{\prime}}}]$ crosses the negative real axis in the complex plane, its square root would experience a sign-reversing branch-cut discontinuity.

To examine the problem more closely, we first consider a single vibrational mode, in which case the matrix expressions reduce back to simple scalar ones. Eq. 12 simplifies to

where $a=(\omega^{2}+\sigma^{2}\Omega^{2})/2\omega\sigma\Omega$. For a positive-curvature, real-frequency mode, $\sigma=1$ and $a\geq 1$, so the expression in the brackets orbits an elliptical trajectory in the complex plane with frequency $\Omega$. At values of $\Omega t=\pi,3\pi,5\pi,\ldots$, the trajectory crosses the negative real axis and the principal square-root $\sqrt{e^{A^{\prime}}}$ changes sign discontinuously. The problem, of course, is related to the leading $e^{i\Omega t}$ like behavior, for which it is clear that $e^{i\Omega t/2}\neq\sqrt{e^{i\Omega t}}$ for some $t$. This suggests that the solution is to factor out the global $e^{i\Omega t}$ dependence of the complex orbit as

We now take the square-root of each factor separately,

where all $\sqrt{\cdots}$ are principal square roots. Because $a\geq 1$ (for real-frequency modes), the quantity $(1-a)/(1+a)$ has an absolute value less than unity and the argument of the right-most square root never crosses the negative real axis. (The quantity in the parentheses makes a simple circular orbit centered about 1.) Its principal square root thus exhibits no branch-cut discontinuity. The global phase of the orbit is carried by the pure exponential factor $e^{-i\sigma\Omega t/2}$, which can be calculated directly because $\sigma$ and $\Omega$ are known. The branch-cut free expression, Eq. 17, also holds for a negative-curvature, imaginary-frequency mode, which has $\sigma=-i$. In this case, $a$ is a purely imaginary number, so that $(1-a)/(1+a)$ is a complex number of unit modulus, and the argument of the square root again does not cross the negative real axis. (This is assuming $t\geq 0$. For $t<0$, one need only calculate $C_{0}(|t|)$ and recognize that $C_{0}(-t)=C_{0}(t)^{*}$.)

This simple factorization trick must now be generalized to the multidimensional case. Starting from Eq. 12, we factor the matrix inverse as

and then calculate its determinant. Recognizing that $\det(\mathbf{A}\mathbf{B})=\det(\mathbf{A})\det(\mathbf{B})$ and $\det(\mathbf{A})=\det(\mathbf{A}^{T})$, we have

These three factors correspond one-to-one with those of Eq. 16, and we need a branch-cut-free square root of each. As in the one-dimensional case, the first factor, $\det(\mathbf{\Lambda}_{+}/2)^{-2}$, is simply a constant and poses no problems. The second factor, $e^{-it\text{Tr}[\mathbf{\Sigma\Omega}]}=\prod_{k}e^{-it\sigma_{k}\Omega_{k}}$, is the product of the global phase factor for each normal mode. For real-frequency modes, this will be a periodic, complex phase. For imaginary-frequency modes, this will be an exponential damping term. As all $\sigma_{k}$ and $\Omega_{k}$ are known, $e^{-it\text{Tr}[\mathbf{\Sigma\Omega}]/2}$ can be calculated directly. The third factor requires more care. At this point in the one-dimensional case, we could directly proceed to the principal square root. However, in the multidimensional case, the determinant in the third factor can still cross the negative real axis, and in general does. A further factorization is necessary, and given that this is a determinant, it is natural to consider the spectrum of the matrix argument. Let $\lambda_{k}$ ($k=1,\ldots,n$) be the eigenvalues of $(1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T})$. Its determinant is then $\lambda_{1}\lambda_{2}\cdots\lambda_{n}$. We assert that each eigenvalue individually does not cross the negative real axis. Therefore, a branch-cut-free square root of the determinant can be constructed by multiplying the principal square root of each eigenvalue respectively. The final determinant is then

In summary, the numerically stable, discontinuity-free vacuum auto-correlation function $C_{0}(t)$ is given by Eqs. 6, 13, 15, and 18.

The general finite-temperature correlation function replaces the vacuum expectation value with a thermal-average trace,

where $\beta=1/kT$ and $Z_{0}(\beta)=\text{Tr}[e^{-\beta H_{0}}]$, the ground state partition function. The finite-temperature spectrum is related to the correlation function by the same Fourier transform as for the zero-temperature case, Eq. 3, with $C_{0}(t)$ replaced by $C(t,\beta)$. In the zero-temperature case, we found a form of the exponential operator in which the vacuum expectation value was simple to evaluate (the disentangled form) and then used the Lie algebra matrix representation to express the exponential in that form. The finite-temperature case proceeds in a parallel fashion. We first find a form in which the trace is simple to evaluate (the diagonalized form, considered below), and then use the same matrix representation to express the exponential product in Eq. 19 in that form.

Consider first the trace of a purely quadratic exponential operator,

where we have assumed that $B_{ij}=C_{ij}$ and $A_{ij}$ is symmetric. The trace is invariant to similarity transformations. Thus, a series of rotation and scaling transformations can diagonalize this operator to

where $\omega_{i}^{\prime 2}$ are the eigenvalues of $\mathbf{G}\mathbf{F}$, for $\mathbf{G}=\mathbf{A}-2\mathbf{B}$ and $\mathbf{F}=\mathbf{A}+2\mathbf{B}$. The parameters of the similarity transformations themselves are not needed. The square of this trace is

which is invariant to the branch choice of $\omega^{\prime}_{i}$. If we let $\mu$ denote the $2n\times 2n$ matrix,

then it can be shown that the squared trace can also be expressed as

Including gradient terms and a constant energy offset adds a contribution of the form $g_{i}(a_{i}+a_{i}^{\dagger})+h$ to the exponential argument. The gradient terms can be eliminated using displacement similarity transformations (which again leave the trace invariant), modifying the total constant term to

Thus, the squared trace of a general symmetric quadratic operator is

The exponential product in Eq. 19 must now be recast into a single exponential operator, which can then be “diagonalized” in the sense above. The trace is first symmetrized as

where $\tau=-it+\beta$. From here, the exponential product is carried out in the finite matrix representation in a similar fashion as the zero-temperature case. The full derivation is shown in Appendix B, which includes the somewhat lengthy expressions for the final branch-cut free and numerically stable result for $C(t,\beta)$. The result is also well-behaved in the $T\rightarrow 0$ ($\beta\rightarrow\infty$) limit.

The anharmonicities of the initial and final state potential energy surfaces both affect the auto-correlation function, the former by modifying the initial wavefunctions and partition function and the latter by modifying the propagation. Considering only the zero-temperature limit, these can be expanded in a perturbation series,

where $\lambda$ is a dummy order-sorting parameter and $V^{\prime}$ is the anharmonic contribution to the upper state surface. The anharmonic term of the initial state surface, $V$, gives corrections to the ground state $|\Psi\rangle$.

Formally, the complete first-order correction to the correlation function is

The first-order initial state term $|\Psi^{(1)}\rangle$ can be computed straightforwardly using time-independent vibrational perturbation theory, expanded in terms of the zeroth-order harmonic oscillator states,

where $\mathbf{m}=\left\{m_{1},m_{2},\ldots,m_{n}\right\}$ are the $n$ quantum numbers of the zeroth-order states. The harmonic cross-correlation functions

can in turn be computed in terms of $C_{0}(t)$ using the multidimensional recurrence relations derived by Fernández and Tipping [4, 27]. On the other hand, the final-state propagator correction $u^{(1)}(t)$ requires a significantly more complex time-dependent treatment, including the renormalization of secular terms [28], which is beyond the scope of this paper. As will be seen below, however, the time-independent initial state term is in many cases the dominant first-order contribution, and the time-dependent final state term can be neglected.²²footnotetext: These recurrence relations exhibit similar numerical instabilities as $C_{0}(t)$ for imaginary-frequency modes, and the same factorizations as those used above provide for stable finite-precision evaluation. We note that the recurrence relations, which provide the ratio of $\langle\mathbf{m}|e^{-iHt}|\mathbf{n}\rangle$ with respect to $C_{0}(t)$, can be used to derive a simple differential equation for $C_{0}(t)$ itself. This differential equation has the feature of having no branch-cut discontinuities. To solve for $C_{0}(t)$ via numerical integration, however, a sufficiently small integration time-step must be used, leading to the same practical circumstance as taking small time-steps to manually correct sign changes in the direct approach.

We anticipate this scenario to occur when there is a large upper state gradient at the vertical geometry. In this case, the leading anharmonic corrections are due to slight changes in the average position $\langle\mathbf{q}\rangle$ of the initial wavefunction $|\Psi\rangle$, which are amplified by the upper state gradient and shift the effective vertical energy. Accounting for cubic anharmoncity, the first-order displacement in each normal coordinate can be evaluated using standard harmonic oscillator matrix elements as

where $V=\frac{1}{6}\sum_{ijk}\phi_{ijk}q_{i}q_{j}q_{k}$ are the ground-state anharmonic terms. The approximate spectral shift is

where $\mathbf{G}$ is the upper-state gradient at the vertical geometry.

Accounting for the vibrational dependence of the electronic transition dipole moment requires consideration of the full dipole auto-correlation function, which for the zero-temperature case is

where we assume the transition dipole function $\mu(\mathbf{q})$ is real. The Franck-Condon approximation truncates the dipole function to its value at the initial state equilibrium geometry $\mu_{e}=\mu(\mathbf{q}=0)$, so that $D_{0}(t)\approx\mu_{e}^{2}C_{0}(t)$. For electronic transitions where $\mu_{e}=0$ by symmetry or $\mu(\mathbf{q})$ varies strongly near $\mathbf{q}=0$, the linear and possibly higher-order $\mathbf{q}$-dependence of $\mu$ (the so-called Herzberg-Teller effect) needs to be accounted for.

Here, we consider just linear dependence,

where $\partial_{i}\mu=\partial\mu/\partial q_{i}|_{e}$ are the equilibrium transition dipole derivatives. The wavefunction $\mu(\mathbf{q})|0\rangle$ can be expanded in terms of the vacuum and one-quantum excited states. The cross-correlation functions of each of these terms are evaluated using the same recurrence relations as in Section II.3. A similar procedure can be used to efficiently combine linear Herzberg-Teller effects with the anharmonic ground state that includes first-order cubic corrections $|\Psi^{(1)}\rangle$.

As with the cubic anharmonicity corrections, the linear dipole dependence will be particularly important when there is a large upper state gradient at the vertical geometry. We can estimate the shift to the effective vertical energy by evaluating the average displacement of $\mu|0\rangle$, which to first order is

i.e., the fractional dipole derivative. The expected spectral energy shift analogous to Eq. 22 is

The absorption spectrum of the dissociative $A$–$X$ transition of F₂ provides a simple one-dimensional example to compare the harmonic correlation and its low-order corrections to the corresponding numerically exact results on a given potential energy curve. We calculated the potential energy and transition dipole curves on a grid of 16 points spaced evenly over $r\in[1.2,1.8]$ Å at the frozen-core CCSDT level of theory [30] with an aug-cc-pVQZ basis set [31] for the $X$ state and its equation-of-motion (EOM) variant [32] for the excited $A$ state. (All electronic structure calculations were performed with the Cfour package [33].) The transition dipole curve was calculated from the geometric mean of the left and right EOM expectation values. Each curve was fitted to an 8${}^{\text{th}}$-order polynomial in the scaled coordinate $x=e^{-r/a}$, with $a=1.0$ Å, from which derivatives were calculated. The potential energy and transition dipole curves are shown in Fig. 1A.

The effects of anharmonicity in the $X$ state surface are illustrated in the top panel of Fig. 1B, where the numerically exact FC spectrum is compared against the harmonic simulation and that which includes first-order ground-state cubic corrections. The first-order spectrum agrees closely with the exact result. The shift of the spectral peak is about $-580$ cm^-1, while the simple cubic-gradient approximation (Eq. 22) predicts $-780$ cm^-1. The corresponding shifts from the linear Herzberg-Teller corrections are shown in the bottom panel of Fig. 1B, where the exact $\mu(r)$-dependent simulation is compared with the harmonic FC simulation and that which includes both first-order cubic and linear Herzberg-Teller corrections. The shift between the simulations with and without Herzberg-Teller corrections is $-360$ cm^-1. The corresponding estimate of this shift (Eq. 24) is $-350$ cm^-1.

The experimentally measured 298 K absorption cross sections [34] are compared to the first-order cubic/linear Herzberg-Teller simulation in Fig. 1C. The absolute absorption cross section is related to the dipole spectral density $D(\omega)$ as

where $g_{e}$ is the electronic degeneracy and $A=2\pi^{2}/3\epsilon_{0}hc\approx 2.6891\times 10^{-18}$ cm${}^{2}/(e\,a_{0})^{2}$. The $A$ state of F₂ has $\Pi$ symmetry, so that $g_{e}=2$. The agreement between the calculated spectrum, which has not been scaled or shifted, and the experimental data is excellent.

As an example of a dissociative system with bound spectator modes, we examined the absorption spectrum of the $\tilde{A}$ and $\tilde{B}$ states of HOCl, both of which dissociate to OH + Cl. The force constants and transition dipole derivatives were determined by a polynomial fit to a grid of 120 ab initio points near the $\tilde{X}$ state equilibrium geometry calculated at the same frozen-core (EOM-)CCSDT/aug-cc-pVQZ level of theory employed for F₂. The experimental absorption cross sections [35, 36, 37] are compared to the simulation with first-order cubic and linear Herzberg-Teller corrections in Fig. 2.

The single largest discrepancy is the total intensity of the $\tilde{B}$–$\tilde{X}$ band. That is, the error arises primarily from the ab initio transition dipole. The accurate shape and breadth of both band systems indicate that the approximate correlation functions and force-fields are accurate. The low-energy shoulder in the experimental spectrum is due to weak absorption from the lowest-lying triplet state [37] and a full treatment of the absolute intensities of this system requires that spin-orbit interactions with the triplet manifold be taken into account [38], well outside the scope of the present study. We also note that finite-temperature effects, calculated at the harmonic approximation, are minor for these systems at 298 K.

Methyleneimine, CH₂NH, is the simplest imine. It has a planar ground-state equilibrium structure [39], and as in the isoelectronic case of ethylene, relaxes to a non-planar geometry with a dihedral angle of 90^∘ between the CNH and CH₂ planes in its first singlet valence excited state [40, 41]. The two equivalent non-planar $C_{s}$ minima of the excited state are connected by a planar $C_{2v}$ transition state with a linear CNH bond angle [40]. Thus, at the ground-state equilibrium geometry, the excited state is best described as a displaced transition state, possessing a large gradient along the CNH bond angle (initially relaxing towards 180^∘) and imaginary frequencies for both the NH torsion and CH₂ out-of-plane angles. Upon $\tilde{A}$–$\tilde{X}$ excitation, the ground state vibrational wavepacket will quickly relax along these degrees of freedom leading to a broad absorption envelope [42].

To simulate the absorption spectrum, we computed the force fields and transition dipole of the $\tilde{X}$ and $\tilde{A}$ states at the frozen-core EOM-CCSD(T)(a)*/aug-cc-pVQZ level of theory [43]. The harmonic and first-order cubic/linear Herzberg-Teller-corrected simulations are compared to the room-temperature experimental spectrum [42] in Fig. 3. The simulation is slightly too high in energy and too low in intensity. Given the relatively small effect of the cubic and Herzberg-Teller corrections, we speculate that this error arises primarily from the vertical electronic energy and equilibrium transition dipole. Indeed, modest shifts of $-450$ cm^-1 to the vertical energy and $+7\%$ to the transition dipole bring the simulation in close agreement with the measurements.

NO₂ exhibits a large change in geometry upon transition to its first excited electronic state, with an equilibrium $\angle$ONO bond angle of 134^∘ in the $\tilde{X}$ state and 102^∘ in the $\tilde{A}$ state [44], leading to an extended Franck-Condon envelope. We simulated the $\tilde{A}-\tilde{X}$ absorption spectrum using EOMIP-CCSD(T)(a)*/aug-cc-pVQZ force fields based on an NO₂ anion reference wavefunction. We assume a vertical transition dipole of $0.302$ a.u. based on the high-level MRCI calculations of Ndengué et al. [44] and supplement this with fractional dipole derivatives calculated at the EOMEE-CCSD/aug-cc-pVQZ level using the neutral $\tilde{X}$ state reference wavefunction. Figure 4A compares the experimental spectrum obtained at 203 K [45, 46] with the harmonic and ground-state-cubic/linear Herzberg-Teller spectra, simulated with a Gaussian FWHM resolution of 500 cm^-1. As with previous examples, the corrections have a significant effect on the harmonic spectrum. The pseudo-Jahn-Teller $\tilde{B}$ and $\tilde{C}$ states, which are neglected here, also make a small contribution to the absorption cross section at the low energy part of the band. While accounting for these states requires a much more elaborate treatment [44], it appears that the cubic/linear Herzberg-Teller simulation of the $\tilde{A}-\tilde{X}$ system alone accounts quantitatively for the gross features of the spectrum.

Because the $\tilde{A}$ state is still being treated harmonically, the auto-correlation function exhibits a regular series of revivals (Fig. 4B). The corresponding high-resolution spectrum (85 cm^-1 FWHM, Fig. 4C) contains extended vibrational progressions. As discussed further below, these high-resolution features are not expected to be quantitatively meaningful at this level of approximation, but have been included here for illustrative purposes.