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

P. Bryan Changala  Center for Astrophysics $|$ Harvard & Smithsonian,
Cambridge, MA, USA<EMAIL_ADDRESS>Nadav Genossar  Joshua H.
Baraban  Department of Chemistry, Ben-Gurion University of the Negev, Beer
Sheva, Israel

###### Abstract

Molecular electronic spectra can be represented in the time domain as auto-
correlation functions of the initial vibrational wavepacket. We present a
derivation of the harmonic vibrational auto-correlation function that is valid
for both real and imaginary harmonic frequencies. The derivation rests on Lie
algebra techniques that map otherwise complicated exponential operator
arithmetic to simpler matrix formulae. The expressions for the zero- and
finite-temperature harmonic auto-correlation functions have been carefully
structured both to be free of branch-cut discontinuities and to remain
numerically stable with finite-precision arithmetic. Simple extensions correct
the harmonic Franck-Condon approximation for the lowest-order anharmonic and
Herzberg-Teller effects. Quantitative simulations are shown for several
examples, including the electronic absorption spectra of F2, HOCl, CH2NH, and
NO2.

##  I Introduction

The vibronic structure of optical, photoionization, and photoelectron spectra
encode the evolution that nuclear wavepackets undergo following electronic
excitation or detachment. Interpreting the dynamics revealed by such spectra
is often aided by first-principles simulations based on two key dynamical
simplifications, the Born-Oppenheimer and Franck-Condon (FC) approximations.
Together these imply that the transition dipole between vibrational levels
belonging to two different electronic states is proportional to the overlap
integral between their respective vibrational wavefunctions. If it is
sufficiently accurate to further approximate the potential energy surfaces as
quadratic expansions about their respective local minima, then the vibrational
structure is that of two sets of (mutually displaced, rotated, and scaled)
harmonic oscillators, the overlap integrals of which are well known [1, 2, 3,
4]. As the only information required in this commonly used framework is the
equilibrium geometries and quadratic force constants of the initial and final
surfaces, there exist a number of popular software packages for computing
harmonic FC spectra using this frequency-domain, state-by-state approach [5,
6, 7, 8]. For large molecules, the number of discrete vibrational states that
contribute to the spectrum grows rapidly, but intelligent selection algorithms
can address this and remain efficient [9, 10, 6].

The frequency-domain approach breaks down when the vertical geometry is highly
displaced from a local minimum (if one exists at all) on the final potential
energy surface. A quadratic expansion centered at the minimum would then be a
poor approximation at the vertical geometry. Although an expansion centered at
the vertical geometry would be more accurate, the large displacement from the
extrapolated effective minimum would still pose a challenge to frequency-
domain methods because of the high vibrational state density at the energy of
the vertical geometry. More fundamentally, the Hessian at the vertical
geometry is not guaranteed to be positive definite: it may possess negative-
curvature (i.e. imaginary-frequency) modes, which generate a continuum of
states in the quadratic approximation. Although one might construct a
frequency-domain method that treats such continua directly (or indirectly
through non-Hermitian techniques), we adopt a different approach.

The problems of large displacements and negative-curvature modes are naturally
addressed by moving to a time-domain picture. Here the FC spectrum is obtained
by the Fourier transform of the auto-correlation function of the initial
vibrational wavepacket evolving on the final-state potential energy surface.
Sophisticated time-domain algorithms, such as the multi-configurational time-
dependent Hartree (MCTDH) method [11] and Chebyshev propagation [12], can
provide highly accurate spectroscopic predictions. These methods, however, are
often computationally expensive, may require special forms of the molecular
Hamiltonian and potential energy surfaces, and typically handle unbound modes
by imposing outgoing boundary conditions using manually tuned complex
absorbing potentials [13, 14, 15]. For spectra that display a broad,
continuous FC envelope, only the short-time, and therefore local, dynamics are
relevant, motivating a time-domain treatment based on perturbation theory, in
the same spirit of standard second-order vibrational perturbation theory
(VPT2) for small-amplitude frequency-domain problems [16]. The zeroth-order
picture of such an approach would be the harmonic auto-correlation function
based on a quadratic expansion of the initial and final potential energy
surfaces (both centered about the initial state equilibrium or vertical
geometry). Analytic expressions for the harmonic correlation function are
known and indeed have been applied to the general vibrational problem before
[17, 18, 19, 20]. These formal expressions, however, suffer from two important
technical problems.

The first problem is branch-cut discontinuities arising from the periodic
motion of bound, real-frequency modes. This analytical issue appears in many
manifestations of the harmonic oscillator correlation function or propagator
(and in the dynamics of oscillators more generally [21, 22]), though it is
often overlooked. This may be because there is a simple numerical fix to the
problem, which has been adopted in prior applications [20, 23, 24].
Nonetheless, it is of interest to derive an analytical solution to the branch-
cut issue for the general, multi-dimensional case, which we present here.
11footnotetext: The procedure is simply to take small enough time-steps to
reliably detect a numerical discontinuity in the correlation function.

The second problem involves unstable finite-precision arithmetic arising from
the exponentially diverging motion of unbound, imaginary-frequency modes.
Although the time-domain expressions used in prior applications are formally
correct for imaginary-frequency modes by analytic continuation, the details of
this procedure require care. Expressions that are well-behaved numerically for
bound, periodic modes may become exponentially divergent for imaginary
frequencies. Unitary evolution guarantees that the final correlation function
is bounded, but this comes about by cancellation of terms that individually
diverge exponentially in time. In any practical numerical implementation,
guaranteeing these cancellations using finite-precision arithmetic requires a
detailed understanding of the analytical structure of the auto-correlation
function. Addressing these two issues is a major focus of this paper.

We first present a derivation of the harmonic correlation function using a
Lie-algebraic approach, for both zero and finite temperatures, which enables
the clear identification of how both the branch-cut and unstable arithmetic
problems arise and how to solve them. Simple extensions beyond the harmonic FC
approximation, including the lowest order anharmonic and Herzberg-Teller
corrections, are then discussed. We demonstrate these results with several
illustrative examples, highlighting applications to quantitative simulations
of the electronic absorption spectra of dissociative surfaces (fluorine, F2,
and hypochlorous acid, HOCl), transition state regions (methyleneimine,
CH2NH), and highly displaced geometries (nitrogen dioxide, NO2).

##  II Theory

The harmonic initial state Hamiltonian, with $\hbar=1$, is

$\displaystyle H_{0}=\sum_{i}^{n}\frac{1}{2}\omega_{i}(p_{i}^{2}+q_{i}^{2}),$
(1)

where $q_{i}$ is the reduced dimensionless normal coordinate of the
$i^{\mathrm{th}}$ vibrational mode ($i=1,\ldots,n$), $p_{i}$ is its conjugate
momentum, and $\omega_{i}$ is its harmonic frequency. The dimensionless
coordinates and momenta satisfy the canonical commutation relation
$[q_{j},p_{k}]=i\delta_{jk}$. In these coordinates, the final state
Hamiltonian expanded to second order about the vertical geometry ($q_{i}=0$)
is

$\displaystyle H$
$\displaystyle=\sum_{i}\frac{1}{2}\omega_{i}p_{i}^{2}+\frac{1}{2}\sum_{ij}K_{ij}q_{i}q_{j}+\sum_{i}G_{i}q_{i}$
$\displaystyle=\frac{1}{2}\mathbf{p}^{T}\mathbf{W}\mathbf{p}+\frac{1}{2}\mathbf{q}^{T}\mathbf{K}\mathbf{q}+\mathbf{G}^{T}\mathbf{q}+V_{0},$
(2)

where $\mathbf{q}$ and $\mathbf{p}$ are $n\times 1$ column vectors containing
the $q_{i}$ and $p_{i}$ operators; $\mathbf{W}$ is a diagonal matrix
containing $\omega_{i}$ along the diagonal; and $V_{0}$, $\mathbf{G}$, and
$\mathbf{K}$ are the final state energy, gradient, and Hessian, respectively,
evaluated at the vertical geometry. The final state Hamiltonian can be written
in terms of its own normal coordinates $\mathbf{Q}$ and momenta $\mathbf{P}$
as

$\displaystyle H$
$\displaystyle=\frac{1}{2}(\mathbf{P}^{T}\mathbf{\Omega}\mathbf{P}+\mathbf{Q}^{T}\mathbf{\Sigma}^{2}\mathbf{\Omega}\mathbf{Q})-\frac{1}{2}\mathbf{G}^{T}\mathbf{K}^{-1}\mathbf{G}+V_{0},$

where

$\displaystyle\mathbf{P}$ $\displaystyle=\mathbf{R}^{-1}\mathbf{p},$
$\displaystyle\mathbf{Q}$
$\displaystyle=\mathbf{R}^{T}(\mathbf{q}-\mathbf{d}),$
$\displaystyle\mathbf{d}$ $\displaystyle=-\mathbf{K}^{-1}\mathbf{G},$
$\displaystyle\mathbf{R}$
$\displaystyle=\mathbf{W}^{-1/2}\mathbf{L}\mathbf{\Omega}^{1/2}.$

The matrix $\mathbf{L}$ is orthogonal and contains the eigenvectors of the
mass-weighted Hessian, i.e.

$\displaystyle(\mathbf{W}^{1/2}\mathbf{K}\mathbf{W}^{1/2})\mathbf{L}=\mathbf{L}\mathbf{Z}$

where $\mathbf{Z}$ contains the real eigenvalues $z_{k}$ along its diagonal.
These may be positive or negative and equal the square of the final-state
harmonic frequencies, which we take to be either positive real or negative
imaginary. That is, let $z_{k}^{1/2}\equiv\sigma_{k}\Omega_{k}$, where
$\Omega_{k}=|z_{k}|^{1/2}$ and $\sigma_{k}=1$ if $z_{k}$ is positive and
$\sigma_{k}=-i$ is $z_{k}$ is negative. The matrices $\mathbf{\Sigma}$ and
$\mathbf{\Omega}$ are diagonal, containing $\sigma_{k}$ and $\Omega_{k}$,
respectively. The matrices $\mathbf{\Lambda}_{\pm}$, which relate the
canonical creation and annihilation operators of the initial and final states,
will appear often below. They are defined as

$\displaystyle\mathbf{\Lambda}_{\pm}$
$\displaystyle=\mathbf{W}^{-1/2}\mathbf{L}\mathbf{\Sigma}^{1/2}\mathbf{\Omega}^{1/2}\pm\mathbf{W}^{1/2}\mathbf{L}\mathbf{\Sigma}^{-1/2}\mathbf{\Omega}^{-1/2}.$

###  II.1 The zero-temperature correlation function

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

$\displaystyle C_{0}(t)=\langle 0|e^{-iHt}|0\rangle,$

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,

$\displaystyle
S_{0}(\omega)=\frac{1}{2\pi}\int_{-\infty}^{\infty}dt\,e^{i(\omega+E_{0})t}C_{0}(t),$
(3)

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

$\displaystyle-iHt=Q\equiv
A_{ij}(a^{\dagger}_{i}a_{j}+\frac{1}{2}\delta_{ij})+B_{ij}a^{\dagger}_{i}a^{\dagger}_{j}+C_{ij}a_{i}a_{j}+f_{i}a_{i}^{\dagger}+g_{i}a_{i}+h,$
(4)

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,

$\displaystyle[Q,Q^{\prime}]=Q^{\prime\prime},$

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

$\displaystyle e^{Q}$
$\displaystyle=e^{A_{ij}(a^{\dagger}_{i}a_{j}+\frac{1}{2}\delta_{ij})+B_{ij}a^{\dagger}_{i}a^{\dagger}_{j}+C_{ij}a_{i}a_{j}+f_{i}a_{i}^{\dagger}+g_{i}a_{i}+h}$
$\displaystyle=e^{h^{\prime}}e^{f^{\prime}_{i}a_{i}^{\dagger}}e^{B^{\prime}_{ij}a_{i}^{\dagger}a_{j}^{\dagger}}e^{A^{\prime}_{ij}(a_{i}^{\dagger}a_{j}+\delta_{ij}/2)}e^{C^{\prime}_{ij}a_{i}a_{j}}e^{g^{\prime}_{i}a_{i}}.$
(5)

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,

$\displaystyle C_{0}(t)$ $\displaystyle=\langle 0|e^{-iHt}|0\rangle$
$\displaystyle=\langle
0|e^{h^{\prime}}e^{f^{\prime}_{i}a_{i}^{\dagger}}e^{B^{\prime}_{ij}a_{i}^{\dagger}a_{j}^{\dagger}}e^{A^{\prime}_{ij}(a_{i}^{\dagger}a_{j}+\delta_{ij}/2)}e^{C^{\prime}_{ij}a_{i}a_{j}}e^{g^{\prime}_{i}a_{i}}|0\rangle$
$\displaystyle=e^{h^{\prime}}\langle
0|e^{A^{\prime}_{ij}(a_{i}^{\dagger}a_{j}+\delta_{ij}/2)}|0\rangle$
$\displaystyle=e^{h^{\prime}}e^{\text{Tr}[\mathbf{A}^{\prime}/2]}$
$\displaystyle=e^{h^{\prime}}\det[e^{\mathbf{A^{\prime}}/2}],$ (6)

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

$\displaystyle M$
$\displaystyle=\left[\begin{array}[]{cccc}0&\mathbf{g}^{T}&\mathbf{f}^{T}&-2h\\\
0&\mathbf{A}&2\mathbf{B}&-\mathbf{f}\\\
0&-2\mathbf{C}&-\mathbf{A}^{T}&\mathbf{g}\\\ 0&0&0&0\end{array}\right].$ (11)

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

$\displaystyle e^{\mathbf{A}^{\prime}}$
$\displaystyle=4\left(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda_{+}}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}\right)^{-1}$
(12)

and

$\displaystyle h^{\prime}$
$\displaystyle=-itV_{0}+\frac{t^{2}}{16}\mathbf{G}^{T}(\mathbf{\Lambda}_{+}-\mathbf{\Lambda}_{-})\mathbf{\Gamma}(\mathbf{\Lambda}_{+}-\mathbf{\Lambda}_{-})^{T}\mathbf{G},$
(13)

where

$\displaystyle\mathbf{\Gamma}$
$\displaystyle=(\zeta^{+}-\zeta^{-})-(\mathbf{\Lambda}_{+}\eta^{+}+\mathbf{\Lambda}_{-}\eta^{-})^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T})^{-1}(\mathbf{\Lambda}_{+}\eta^{+}+\mathbf{\Lambda}_{-}\eta^{-}),$
(14)

and

$\displaystyle e^{\pm}$ $\displaystyle=e^{\pm
it\mathbf{\Sigma}\mathbf{\Omega}},$ $\displaystyle\eta^{\pm}$
$\displaystyle=\frac{e^{\pm it\mathbf{\Sigma}\mathbf{\Omega}}-1}{\pm
it\mathbf{\Sigma}\mathbf{\Omega}},$ $\displaystyle\zeta^{\pm}$
$\displaystyle=\frac{e^{\pm it\mathbf{\Sigma}\mathbf{\Omega}}\mp
it\mathbf{\Sigma}\mathbf{\Omega}-1}{-t^{2}(\mathbf{\Sigma\Omega})^{2}}.$

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,

$\displaystyle\mathbf{\Gamma}$
$\displaystyle\rightarrow\zeta^{+}-\eta^{+}\mathbf{\Lambda_{+}}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T})^{-1}\mathbf{\Lambda}_{+}\eta^{+}$
$\displaystyle=\zeta^{+}-\eta^{+}e^{-}\eta^{+}$
$\displaystyle=\zeta^{+}-\eta^{+}\eta^{-}$
$\displaystyle=\zeta^{+}-(\zeta^{+}+\zeta^{-})$ $\displaystyle=-\zeta^{-}.$

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,

$\displaystyle\left(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda_{+}}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}\right)^{-1}$
$\displaystyle=\left(1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}\right)^{-1}\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}$
$\displaystyle=\mathbf{\Lambda}^{\prime}_{+}\left(1-e^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}\mathbf{\Lambda}^{\prime}_{+}\right)^{-1}e^{-}\mathbf{\Lambda}_{+}^{-1}$
$\displaystyle=\mathbf{\Lambda}^{\prime}_{+}e^{-}\left(1-\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}\mathbf{\Lambda}^{\prime}_{+}e^{-}\right)^{-1}\mathbf{\Lambda}_{+}^{-1}$
$\displaystyle=\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}\left(1-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}\right)^{-1},$

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

$\displaystyle\mathbf{\Gamma}$
$\displaystyle=-\eta^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}\mathbf{\Lambda}^{\prime}_{+}(1-e^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}\mathbf{\Lambda}^{\prime}_{+})^{-1}\eta^{-}$
$\displaystyle\qquad-\left[\eta^{-}\mathbf{\Lambda}_{-}^{T}(1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T})^{-1}\mathbf{\Lambda}^{\prime}_{+}\eta^{-}+\text{transpose}\right]$
$\displaystyle\qquad-2\zeta^{-}-\eta^{-}\mathbf{\Lambda}_{-}^{T}(1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T})^{-1}\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}\eta^{-}.$
(15)

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

$\displaystyle e^{A^{\prime}}$ $\displaystyle=\left[\cos(\sigma\Omega
t)+ia\sin(\sigma\Omega t)\right]^{-1},$

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

$\displaystyle e^{A^{\prime}}=\frac{2}{1+a}e^{-i\sigma\Omega
t}\left(1+\frac{1-a}{1+a}e^{-2i\sigma\Omega t}\right)^{-1}.$ (16)

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

$\displaystyle e^{A^{\prime}/2}=\sqrt{\frac{2}{1+a}}e^{-i\sigma\Omega
t/2}\sqrt{\left(1+\frac{1-a}{1+a}e^{-2i\sigma\Omega t}\right)^{-1}},$ (17)

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

$\displaystyle e^{\mathbf{A}^{\prime}}$
$\displaystyle=4(1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T})^{-1}\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}$

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

$\displaystyle\det(e^{\mathbf{A}^{\prime}})$
$\displaystyle=\det(\mathbf{\Lambda}_{+}/2)^{-2}e^{-it\text{Tr}[\mathbf{\Sigma\Omega}]}\det(1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T})^{-1}.$

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

$\displaystyle\det(e^{\mathbf{A}^{\prime}/2})=\det(\mathbf{\Lambda}_{+}/2)^{-1}e^{-it\text{Tr}[\mathbf{\Sigma\Omega}]/2}\prod_{k}\sqrt{\lambda_{k}^{-1}}.$
(18)

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

###  II.2 The finite-temperature correlation function

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

$\displaystyle C(t,\beta)$
$\displaystyle=\frac{1}{Z_{0}(\beta)}\text{Tr}\left[e^{(+it-\beta)H_{0}}e^{-itH}\right],$
(19)

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,

$\displaystyle\text{Tr}\big{[}\exp[A_{ij}(a^{\dagger}_{i}a_{j}+\delta_{ij}/2)+B_{ij}(a_{i}a_{j}+a^{\dagger}_{i}a^{\dagger}_{j})]\big{]},$

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

$\displaystyle\text{Tr}[e^{\omega^{\prime}_{i}(a^{\dagger}_{i}a_{i}+1/2)}]$
$\displaystyle=\prod_{i}\frac{e^{\omega^{\prime}_{i}/2}}{1-e^{\omega^{\prime}_{i}}},$

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

$\displaystyle\text{Tr}[\cdots]^{2}$
$\displaystyle=\prod_{i}\frac{1}{2}\frac{1}{\cosh\omega_{i}^{\prime}-1},$

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

$\displaystyle\mu$
$\displaystyle=\left[\begin{array}[]{cc}\mathbf{A}&2\mathbf{B}\\\
-2\mathbf{B}&-\mathbf{A}\end{array}\right],$

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

$\displaystyle\text{Tr}[\cdots]^{2}$
$\displaystyle=(-1)^{n}\det\left[e^{\mu}-\mathbf{1}\right]^{-1}.$

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

$\displaystyle
h^{\prime}=-\mathbf{g}^{T}\left(\mathbf{A}+2\mathbf{B}\right)^{-1}\mathbf{g}+h.$
(20)

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

$\displaystyle\text{Tr}[\cdots]^{2}$
$\displaystyle=(-1)^{n}\det|e^{\mu}-\mathbf{1}|^{-1}\times e^{2h^{\prime}}.$

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

$\displaystyle\text{Tr}\left[e^{(+it-\beta)H_{0}}e^{-itH_{1}}\right]$
$\displaystyle=\text{Tr}\left[e^{-\tau H_{0}}e^{-itH_{1}}\right]$
$\displaystyle=\text{Tr}\left[e^{-\tau H_{0}/2}e^{-itH_{1}}e^{-\tau
H_{0}/2}\right],$ (21)

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.

###  II.3 Simple anharmonic corrections

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,

$\displaystyle|\Psi\rangle$
$\displaystyle=|0\rangle+\lambda|\Psi^{(1)}\rangle+\lambda^{2}|\Psi^{(2)}\rangle\cdots$
$\displaystyle e^{-i(H+V^{\prime})t}$ $\displaystyle=e^{-iHt}(1+\lambda
u^{(1)}(t)+\lambda^{2}u^{(2)}(t)\cdots),$

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

$\displaystyle\lambda C_{0}^{(1)}(t)$ $\displaystyle=\lambda\left(\langle
0|e^{-iHt}|\Psi^{(1)}\rangle+\langle\Psi^{(1)}|e^{-iHt}|0\rangle+\langle
0|e^{-iHt}u^{(1)}(t)|0\rangle\right)$ $\displaystyle=\lambda\langle
0|e^{-iHt}\left(2|\Psi^{(1)}\rangle+u^{(1)}(t)|0\rangle\right).$

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,

$\displaystyle|\Psi^{(1)}\rangle=\sum_{\mathbf{m}}c_{\mathbf{m}}|\mathbf{m}\rangle,$

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

$\displaystyle\langle\mathbf{m}|e^{-iHt}|0\rangle,$

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.22footnotetext: 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

$\displaystyle\langle q_{i}\rangle$
$\displaystyle=-\frac{1}{4\omega_{i}}\sum_{j}\phi_{ijj},$

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

$\displaystyle\Delta E_{\text{cubic}}$
$\displaystyle\approx\mathbf{G}\cdot\langle\mathbf{q}\rangle$
$\displaystyle=-\frac{1}{4}\sum_{ij}\frac{G_{i}\phi_{ijj}}{\omega_{i}},$ (22)

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

###  II.4 Linear Herzberg-Teller corrections

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

$\displaystyle D_{0}(t)$ $\displaystyle=\langle
0|\mu(\mathbf{q})e^{-iHt}\mu(\mathbf{q})|0\rangle,$ (23)

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,

$\displaystyle\mu(\mathbf{q})=\mu_{e}+\sum_{i}\partial_{i}\mu\,q_{i},$

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

$\displaystyle\langle q_{i}\rangle$ $\displaystyle=\frac{\langle 0|\mu
q_{i}\mu|0\rangle}{\langle 0|\mu^{2}|0\rangle}$
$\displaystyle\approx\frac{\partial_{i}\mu}{\mu_{e}},$

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

$\displaystyle\Delta E_{\text{HT}}$
$\displaystyle\approx\mathbf{G}\cdot\langle\mathbf{q}\rangle$
$\displaystyle=\sum_{i}G_{i}\frac{\partial_{i}\mu}{\mu_{e}}.$ (24)

##  III Examples

We now explore various applications of the harmonic correlation function and
its low-order anharmonic and Herzberg-Teller corrections as implemented in the
Nitrogen Python package [29]. Focus is placed on systems for which the time-
domain approach is most appropriate. The correctness of the implementation for
bound harmonic calculations has been verified by comparison with frequency-
domain programs [6] and for unbound problems by comparison with numerical
wavepacket propagation on low-dimensional model potentials.

###  III.1 Dissociative excited states of F2 and HOCl

The absorption spectrum of the dissociative $A$–$X$ transition of F2 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.

Figure 1: The $A$–$X$ band of F2. (A) The potential energy curves (top panel)
and transition dipole (bottom panel) of the $X$ and $A$ states show the ab
initio points (dots) and the fitted curves (solid line). The ground
vibrational wavefunction is plotted over the $X$ state curve. (B) The top
panel compares the exact Franck-Condon simulation (black, solid) with the
harmonic (red, dashed) and first-order cubic (blue, dot-dashed) simulations.
The bottom panel shows the exact $\mu(r)$-dependent simulation (black, solid)
with the harmonic (red, dashed) and first-order cubic plus linear Herzberg-
Teller (blue, dot-dashed) simulations. (C) The calculated first-order cubic
plus linear Herzberg-Teller absolute absorption cross sections (solid) are
compared with the 298 K experimental spectrum [34] (dots).

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

$\displaystyle\sigma(\omega)=A\times g_{e}\times\omega D(\omega),$

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

Figure 2: The absorption spectrum of the $\tilde{A}$–$\tilde{X}$ and
$\tilde{B}$–$\tilde{X}$ bands of HOCl. The EOM-CCSDT/aug-cc-pVQZ simulation
includes first-order cubic and linear Herzberg-Teller corrections. The
contributions from the $\tilde{A}$ (red, dashed) and $\tilde{B}$ (blue, dot-
dashed) states are shown together with their sum (black, solid). The
experimental data points (dots) are at 298 K [35].

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

###  III.2 Accessing isomerization barriers: CH2NH

Methyleneimine, CH2NH, 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 CH2 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 CH2 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.

Figure 3: The absorption spectrum of the $\tilde{A}$–$\tilde{X}$ band of
CH2NH. The room-temperature experimental spectrum [42] is shown with the
harmonic (black, dashed) and first-order cubic/linear Herzberg-Teller-
corrected (black, solid) simulations at the EOM-CCSD(T)(a)*/aug-cc-pVQZ level
of theory. (Simulated resolution = 1050 cm-1 FWHM.) A scaled simulation with
the vertical energy decreased by 450 cm-1 and the transition dipole increased
by 7% is also shown (red, solid).

###  III.3 Highly displaced vertical geometry: NO2

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

Figure 4: The absorption spectrum of the $\tilde{A}$–$\tilde{X}$ band of NO2.
(A) The experimental spectrum [45, 46] at 203 K (thin solid) is shown with the
harmonic (red, dashed) and first-order cubic/linear Herzberg-Teller-corrected
(blue, dot-dashed) simulations (FWHM = 500 cm-1). The absolute value of the
harmonic auto-correlation function shows periodic wavepacket revivals (B),
leading to regular vibrational progressions in a high-resolution (FWHM = 85
cm-1) simulation (C).

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.

##  IV Discussion

The parameters of a vibronic spectrum are naturally partitioned into nearly
independent contributions from the vertical electronic energy, the transition
dipole, and their respective force-fields or derivatives. As the examples
above show, when the potential energy surfaces of interest are approximately
quadratic, the spectral contours can be accurately reproduced by analytical
harmonic correlation functions augmented by low-order anharmonic and Herzberg-
Teller corrections. In this limit, the remaining disagreement with
measurements appears to be dominated by purely electronic errors in the
vertical energy and oscillator strength. We stress that in dissociative or
quasi-unbound spectra, what might normally be categorized separately as
“energy” versus “intensity” effects are intrinsically convoluted. An accurate
prediction of even basic descriptors, such as the peak spectral wavelength,
require a balanced treatment of both. Indeed, this convolution poses a
challenge to extensive quantum chemical benchmarking of vertical energies and
oscillator strengths. A relatively inexpensive vibrational treatment as
presented here may help bridge this gap, at least for those systems that
possess simple enough potential energy surfaces in the vertical region.

The applications above have focused on systems with unstructured spectra, the
envelopes of which are determined by the initial decay of the correlation
function. This is not an inherent limitation of the time-dependent approach.
In specific cases of conventional bound-to-bound transitions where the regions
of the PES accessible to the propagating wavepacket are accurately described
in the harmonic approximation, the long-time behavior of the correlation
function is reliable, and in all cases the frequency-domain and time-domain
results will of course be equivalent. When the vertical geometry is highly
displaced from the upper-state minimum, however, the long-time behavior of the
harmonic correlation function will not be accurate. The $\tilde{A}$ state of
NO2 examined above is a good example of this. Here, the vertical geometry is
displaced from the upper state minimum by approximately $\Delta Q=10$ in terms
of the dimensionless bending normal coordinate. Although the harmonic
correlation function exhibits partial revivals at the vibrational periods of
the harmonic normal modes, these revivals are not meaningful given that
between them the wavepacket traverses an extended region of the PES poorly
described by the harmonic expansion. The same problem effectively occurs in a
harmonic frequency-domain simulation. Accurate long-time simulations may be
possible using time-dependent perturbation theory if the relevant region of
the PES is well approximated with low-order anharmonic corrections. Bound
vibrational modes, in particular, will probably require some sort of
renormalized perturbative expansion of the correlation function [28, 47].
33footnotetext: In fact, the use of a vertical force field, rather than that
of the upper-state minimum, is in some sense a crude form of renormalization.

Finally, while we have focused on electronic absorption spectra, these methods
are applicable to other techniques where FC or low-order Herzberg-Teller
approximations are routinely employed, including photodetachment,
photoionization, and photoelectron spectroscopy [48], as well as non-radiative
processes [49, 50, 51]. It may also be possible to incorporate directly
rotational effects such as origin shifts or axis-switching [52], which are
important mostly for light molecules, especially small hydrides. Other terms
in the rovibrational Hamiltonian, such as Coriolis coupling, should naturally
be considered as higher-order perturbative anharmonic corrections to this
framework are developed.

##  V Conclusions

We have presented a new derivation of the multidimensional harmonic oscillator
auto-correlation function focused on the analytical and numerical issues
associated with branch-cuts and stable finite-precision arithmetic. Together
with the lowest-order anharmonic and Herzberg-Teller corrections, these
results provide a simple means to simulating vibronic spectra that access
unbound or negative curvature potential energy surfaces. The development of
higher-order time-dependent perturbation theory corrections presents a
promising avenue for future work.

###### Acknowledgements.

This project has received funding from the European Research Council (ERC)
under the European Union’s Horizon 2020 research and innovation programme
(grant agreement No. 848668). This work was also supported by the Israel
Science Foundation (ISF), grant No. 194/20, Grant No. 2020724 from the United
States-Israel Binational Science Foundation (BSF) and the US National Science
Foundation (NSF award PHY-2110489).

## Conflict of Interest

The authors declare no conflicts of interest.

## Author’s Contributions

All authors contributed equally to this work.

## Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Derivation of Eqs. 12 and 13

We first consider the purely quadratic sub-algebra represented by the matrices
$\mu$, which have the general form

$\displaystyle\mu$
$\displaystyle=\left(\begin{array}[]{cc}\mathbf{A}&2\mathbf{B}\\\
-2\mathbf{C}&-\mathbf{A}^{T}\end{array}\right).$

This is simply the center block of the full matrix representation, Eq. 11. For
an exponential argument $-iHt$, this simplifies to

$\displaystyle\mu$
$\displaystyle=\frac{-it}{2}\left(\begin{array}[]{c@{\quad}c}\mathbf{W}+\mathbf{K}&-\mathbf{W}+\mathbf{K}\\\
\mathbf{W}-\mathbf{K}&-\mathbf{W}-\mathbf{K}\end{array}\right),$

where $\mathbf{W}$ and $\mathbf{K}$ are defined in Eq. 2. This matrix has
right eigenvectors

$\displaystyle\mathcal{L}$
$\displaystyle=\frac{1}{2}\left(\begin{array}[]{cc}\mathbf{\Lambda}_{+}&-\mathbf{\Lambda}_{-}\\\
-\mathbf{\Lambda}_{-}&\mathbf{\Lambda}_{+}\end{array}\right)$

and left eigenvectors

$\displaystyle\mathcal{L}^{-1}$
$\displaystyle=\frac{1}{2}\left(\begin{array}[]{cc}\mathbf{\Lambda}_{+}^{T}&\mathbf{\Lambda}_{-}^{T}\\\
\mathbf{\Lambda}_{-}^{T}&\mathbf{\Lambda}_{+}^{T}\end{array}\right),$

so that $\mu$ can be diagonalized as

$\displaystyle\mu$
$\displaystyle=-it\times\mathcal{L}\left(\begin{array}[]{cc}\mathbf{\Sigma}\mathbf{\Omega}&0\\\
0&-\mathbf{\Sigma}\mathbf{\Omega}\end{array}\right)\mathcal{L}^{-1}.$

This permits its exponential to be simply expressed as

$\displaystyle e^{\mu}$
$\displaystyle=\mathcal{L}\left(\begin{array}[]{cc}e^{-}&0\\\
0&e^{+}\end{array}\right)\mathcal{L}^{-1},$ (27)

where $e^{\pm}\equiv\exp[\pm i\mathbf{\Sigma\Omega}t]$.

The full matrix representation of the general quadratic Hamiltonian, including
gradient and constant terms, is

$\displaystyle M$
$\displaystyle=\left(\begin{array}[]{ccc}0&\left[\begin{array}[]{cc}\mathbf{g}^{T}&\mathbf{f}^{T}\end{array}\right]&-2h\\\
0&\mu&\left[\begin{array}[]{c}-\mathbf{f}\\\ \mathbf{g}\end{array}\right]\\\
0&0&0\end{array}\right).$

Repeated multiplication shows that higher powers of this matrix are

$\displaystyle M^{k}$
$\displaystyle=\left(\begin{array}[]{c@{\quad}c@{\quad}c}0&\left[\begin{array}[]{cc}\mathbf{g}^{T}&\mathbf{f}^{T}\end{array}\right]\mu^{k-1}&\left[\begin{array}[]{cc}\mathbf{g}^{T}&\mathbf{f}^{T}\end{array}\right]\mu^{k-2}\left[\begin{array}[]{c}-\mathbf{f}\\\
\mathbf{g}\end{array}\right]\\\
0&\mu^{k}&\mu^{k-1}\left[\begin{array}[]{c}-\mathbf{f}\\\
\mathbf{g}\end{array}\right]\\\ 0&0&0\end{array}\right)$

for $k\geq 2$. By direct summation, the matrix exponential is therefore

$\displaystyle e^{M}$
$\displaystyle=\left(\begin{array}[]{c@{\qquad}c@{\qquad}c}1&\left[\begin{array}[]{cc}\mathbf{g}^{T}&\mathbf{f}^{T}\end{array}\right]\mu^{-1}(e^{\mu}-1)&\left[\begin{array}[]{cc}\mathbf{g}^{T}&\mathbf{f}^{T}\end{array}\right]\mu^{-2}(e^{\mu}-\mu-1)\left[\begin{array}[]{c}-\mathbf{f}\\\
\mathbf{g}\end{array}\right]-2h\\\
0&e^{\mu}&\mu^{-1}(e^{\mu}-1)\left[\begin{array}[]{c}-\mathbf{f}\\\
\mathbf{g}\end{array}\right]\\\ 0&0&1\end{array}\right).$ (37)

We now consider the disentangled operator,
$e^{h^{\prime}}e^{f^{\prime}_{i}a_{i}^{\dagger}}e^{B^{\prime}_{ij}a_{i}^{\dagger}a_{j}^{\dagger}}e^{A^{\prime}_{ij}(a_{i}^{\dagger}a_{j}+\delta_{ij}/2)}e^{C^{\prime}_{ij}a_{i}a_{j}}e^{g^{\prime}_{i}a_{i}}$.
The matrix representation of each exponential factor can be evaluated using
Eq. 37. We assume a Hermitian Hamiltonian, for which
$\mathbf{B}^{\prime}=\mathbf{C}^{\prime}$,
$\mathbf{A}^{\prime}=\mathbf{A}^{\prime T}$, and
$\mathbf{g}^{\prime}=\mathbf{f}^{\prime}$. After matrix multiplication, the
final product is

$\displaystyle\left(\begin{array}[]{c@{\qquad}c@{\qquad}c}1&\left[\begin{array}[]{c@{\quad}c}\mathbf{g}^{\prime
T}-2\mathbf{g}^{\prime
T}e^{-\mathbf{A}^{\prime}}\mathbf{B}^{\prime}&\mathbf{g}^{\prime
T}e^{-\mathbf{A}^{\prime}}\end{array}\right]&\mathbf{g}^{\prime
T}e^{-\mathbf{A}^{\prime}}\mathbf{g}^{\prime}-2h^{\prime}\\\
0&\left[\begin{array}[]{c@{\quad}c}e^{\mathbf{A}^{\prime}}-4\mathbf{B}^{\prime}e^{-\mathbf{A}^{\prime}}\mathbf{B}^{\prime}&2\mathbf{B}^{\prime}e^{-\mathbf{A}^{\prime}}\\\
-2e^{-\mathbf{A}^{\prime}}\mathbf{B}^{\prime}&e^{-\mathbf{A}^{\prime}}\end{array}\right]&\left[\begin{array}[]{c}-\mathbf{g}^{\prime}+2\mathbf{B}^{\prime}e^{-\mathbf{A}^{\prime}}\mathbf{g}^{\prime}\\\
e^{-\mathbf{A}^{\prime}}\mathbf{g}^{\prime}\end{array}\right]\\\
0&0&1\end{array}\right).$ (46)

Equating Eq. 37 and Eq. 46 allows us to solve for the disentangled parameters
block-by-block. For example, the lower-right element of the central $2n\times
2n$ block yields

$\displaystyle e^{-\mathbf{A}^{\prime}}$
$\displaystyle=\frac{1}{4}\left(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda_{+}}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}\right),$

where $e^{\mu}$ has been expanded using Eq. 27. The matrix inverse of this
expression gives Eq. 12. Similarly, equating the upper-right elements and
solving for $h^{\prime}$ gives Eq. 13.

##  Appendix B Derivation of the finite-temperature trace

Beginning with the symmetric form of the thermal trace, Eq. 21, the matrix
representation of each factor is directly exponentiated and then multiplied.
The final product is

$\displaystyle\left(\begin{array}[]{c@{\quad}c@{\quad}c}1&\frac{-it}{\sqrt{2}}\left[\begin{array}[]{cc}\mathbf{G}^{T}&\mathbf{G}^{T}\end{array}\right]\mathcal{L}\underline{\eta}\mathcal{L}^{-1}\underline{\xi}^{1/2}&\frac{-t^{2}}{2}\left[\begin{array}[]{cc}\mathbf{G}^{T}&\mathbf{G}^{T}\end{array}\right]\mathcal{L}\underline{\zeta}\mathcal{L}^{-1}\left[\begin{array}[]{c}-\mathbf{G}\\\
\mathbf{G}\end{array}\right]+2itV_{0}\\\
0&\underline{\xi}^{1/2}\mathcal{L}\underline{e}\mathcal{L}^{-1}\underline{\xi}^{1/2}&\frac{-it}{\sqrt{2}}\underline{\xi}^{1/2}\mathcal{L}\underline{\eta}\mathcal{L}^{-1}\left[\begin{array}[]{c}-\mathbf{G}\\\
\mathbf{G}\end{array}\right]\\\ 0&0&1\end{array}\right),$ (56)

where $\xi=\exp[\mathbf{W}\tau]$ and

$\displaystyle\underline{e}$
$\displaystyle=\left[\begin{array}[]{cc}e^{-}&0\\\ 0&e^{+}\end{array}\right],$
$\displaystyle\underline{\eta}$
$\displaystyle=\left[\begin{array}[]{cc}\eta^{-}&0\\\
0&\eta^{+}\end{array}\right],$ $\displaystyle\underline{\zeta}$
$\displaystyle=\left[\begin{array}[]{cc}\zeta^{-}&0\\\
0&\zeta^{+}\end{array}\right],$ $\displaystyle\underline{\xi}^{1/2}$
$\displaystyle=\left[\begin{array}[]{cc}\xi^{-1/2}&0\\\
0&\xi^{1/2}\end{array}\right].$

We now note a useful relation that allows us to evaluate $h^{\prime}$, defined
in Eq. 20, using only elements of the exponentiated representation, Eq. 37,
not the quadratic operator argument itself. First, subtract identity from
$e^{\mu}$ and invert, then sandwich this between the upper and right edges of
the exponential representation. Finally, subtract the upper-right corner
element and divide by 2. The result is

$\displaystyle\frac{1}{2}\left[\begin{array}[]{cc}\mathbf{g}^{T}&\mathbf{g}^{T}\end{array}\right]\mu^{-1}\left[\begin{array}[]{c}-\mathbf{g}\\\
\mathbf{g}\end{array}\right]+h$
$\displaystyle=-\mathbf{g}^{T}(\mathbf{A}+2\mathbf{B})^{-1}\mathbf{g}+h$
$\displaystyle=h^{\prime}.$

Applying this procedure to Eq. 56 yields

$\displaystyle h^{\prime}$
$\displaystyle=-itV_{0}+\frac{t^{2}}{4}\left[\begin{array}[]{cc}\mathbf{G}^{T}&\mathbf{G}^{T}\end{array}\right]\mathcal{L}\left\\{\underline{\zeta}-\underline{\eta}\left(\underline{e}-\mathcal{L}^{-1}\underline{\xi}^{-1}\mathcal{L}\right)^{-1}\underline{\eta}\right\\}\mathcal{L}^{-1}\left[\begin{array}[]{c}-\mathbf{G}\\\
\mathbf{G}\end{array}\right].$

The regularization of the gradient terms starts with expanding this expression
for $h^{\prime}$. The matrix inversion is handled block-wise assuming the
bottom-right block and its respective Schur complement are invertible. Then
the matrix multiplication is carried out and terms are recollected as single-
block expressions. After lengthy algebra, the result has the same form as Eq.
13,

$\displaystyle
h^{\prime}=-itV_{0}+\frac{t^{2}}{16}\mathbf{G}^{T}(\mathbf{\Lambda}_{+}-\mathbf{\Lambda}_{-})\mathbf{\Gamma}(\mathbf{\Lambda}_{+}-\mathbf{\Lambda}_{-})^{T}\mathbf{G},$

where

$\displaystyle\mathbf{\Gamma}$
$\displaystyle=(\zeta^{+}-\zeta^{-})-\theta_{+}^{T}\bar{D}^{-1}\theta_{+}$
(57a)
$\displaystyle\qquad+\left(\theta_{-}-\bar{C}^{T}\bar{D}^{-1}\theta_{+}\right)^{T}\xi^{-1/2}\left(\bar{A}+\xi^{-1/2}\bar{C}^{T}\bar{D}^{-1}\bar{C}\xi^{-1/2}\right)^{-1}\xi^{-1/2}\left(\theta_{-}-\bar{C}^{T}\bar{D}^{-1}\theta_{+}\right)$
(57b)

and

$\displaystyle\theta_{\pm}$
$\displaystyle=\mathbf{\Lambda}_{\pm}\eta^{+}+\mathbf{\Lambda}_{\mp}\eta^{-}$
$\displaystyle\bar{A}$
$\displaystyle=\xi^{-1/2}(\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{-}^{T})\xi^{-1/2}-4$
$\displaystyle\bar{C}$
$\displaystyle=\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{-}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{+}^{T}$
$\displaystyle\bar{D}$
$\displaystyle=\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1}$

In the zero temperature limit, $\xi^{-1}\rightarrow 0$, only the first line of
$\mathbf{\Gamma}$ remains (Eq. 57a), and it reduces to the vacuum expectation
value expression above. Its regularization is handled similarly,

$\displaystyle(\zeta^{+}-\zeta^{-})-\theta_{+}^{T}\bar{D}^{-1}\theta_{+}$
$\displaystyle=(1)+(2)+(3)$ $\displaystyle(1)$
$\displaystyle=\eta^{+}\eta^{-}-\eta^{+}\mathbf{\Lambda}_{+}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{+}\eta^{+}$
$\displaystyle=\eta^{-}e^{+}\left[1-\mathbf{\Lambda}_{+}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{+}e^{+}\right]\eta^{-}$
$\displaystyle=\eta^{-}e^{+}\left[1-(1-e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+})^{-1}\right]\eta^{-}$
$\displaystyle=\eta^{-}e^{+}\left[(-e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+})(1-e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+})^{-1}\right]\eta^{-}$
$\displaystyle=-\eta^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}(1-e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+})^{-1}\eta^{-}$
$\displaystyle(2)$
$\displaystyle=-2\zeta^{-1}-\eta^{-}\mathbf{\Lambda}_{-}^{T}\bar{D}^{-1}\mathbf{\Lambda}_{-}\eta^{-}$
$\displaystyle(3)$
$\displaystyle=-\eta^{+}\mathbf{\Lambda}_{+}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{-}\eta^{-}+\text{transpose}$
$\displaystyle=-\eta^{-}(1-\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}e^{-})^{-1}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}\eta^{-}+\text{transpose}$

The remaining terms are

$\displaystyle\theta_{-}-\bar{C}^{T}\bar{D}^{-1}\theta_{+}$
$\displaystyle=(1)+(2)+(3)+(4)$ $\displaystyle(1)$
$\displaystyle=\mathbf{\Lambda}_{-}\eta^{+}-\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{+}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{+}\eta^{+}$
$\displaystyle=\mathbf{\Lambda}_{-}\left[1-(1-\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}e^{-})^{-1}\right]\eta^{+}$
$\displaystyle=\mathbf{\Lambda}_{-}\left[(1-\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}e^{-})^{-1}(-1)\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}e^{-}\right]\eta^{+}$
$\displaystyle=-\mathbf{\Lambda}_{-}(1-\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}e^{-})^{-1}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}\eta^{-}$
$\displaystyle(2)$
$\displaystyle=\mathbf{\Lambda}_{+}\eta^{-}+\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{-}^{T}\bar{D}^{-1}\mathbf{\Lambda}_{-}\eta^{-}$
$\displaystyle(3)$
$\displaystyle=-\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{+}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{-}\eta^{-}$
$\displaystyle=-\mathbf{\Lambda}_{-}(1-\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}e^{-})^{-1}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}\eta^{-}$
$\displaystyle(4)$
$\displaystyle=\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{-}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{+}\eta^{+}$
$\displaystyle=\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{-}^{T}\mathbf{\Lambda}^{\prime}_{+}(1-e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+})^{-1}\eta^{-}$

and

$\displaystyle\bar{A}+\xi^{-1/2}$
$\displaystyle\bar{C}^{T}\bar{D}^{-1}\bar{C}\xi^{-1/2}=(1)+(2)+(3)$
$\displaystyle(1)$
$\displaystyle=\xi^{-1/2}\left[-\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{-}^{T}+\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{+}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{-}^{T}\right]\xi^{-1/2}$
$\displaystyle=\xi^{-1/2}\mathbf{\Lambda}_{-}e^{+}\left[-1+\mathbf{\Lambda}_{+}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{+}e^{+}\right]\mathbf{\Lambda}_{-}^{T}\xi^{-1/2}$
$\displaystyle=\xi^{-1/2}\mathbf{\Lambda}_{-}e^{+}\left[-1+(1-e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+})^{-1}\right]\mathbf{\Lambda}_{-}^{T}\xi^{-1/2}$
$\displaystyle=\xi^{-1/2}\mathbf{\Lambda}_{-}e^{+}\left[e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}(1-e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+})^{-1}\right]\mathbf{\Lambda}_{-}^{T}\xi^{-1/2}$
$\displaystyle=\xi^{-1/2}\mathbf{\Lambda}_{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}(1-e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+})^{-1}\mathbf{\Lambda}_{-}^{T}\xi^{-1/2}$
$\displaystyle(2)$
$\displaystyle=\xi^{-1/2}\left[\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{+}^{T}+\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{-}^{T}\bar{D}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{+}^{T}\right]\xi^{-1/2}-4$
$\displaystyle(3)$
$\displaystyle=-\xi^{-1/2}\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{+}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{+}^{T}\xi^{-1/2}+\text{transpose}$
$\displaystyle=-\xi^{-1/2}\mathbf{\Lambda}_{-}(1-\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\mathbf{\Lambda}^{\prime}_{+}e^{-})^{-1}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{+}^{T}\xi^{-1/2}+\text{transpose}$

The remaining factor of the thermal correlation function contains the purely
quadratic part of the thermal trace and the initial state partition function,
$Z_{0}$. Together, these are

$\displaystyle(*)$
$\displaystyle\equiv(-1)^{n}\frac{1}{\det[e^{\mu}-\mathbf{1}]Z_{0}^{2}}$
$\displaystyle=\frac{(-1)^{n}}{Z_{0}^{2}}\det\left[\underline{\xi}^{1/2}\mathcal{L}\underline{e}\mathcal{L}^{-1}\underline{\xi}^{1/2}-\mathbf{1}\right]^{-1}$
$\displaystyle=\frac{(-1)^{n}}{Z_{0}^{2}}\det\left[\mathcal{L}\underline{e}\mathcal{L}^{-1}-\underline{\xi}^{-1}\right]^{-1}$

The harmonic partition function is

$\displaystyle Z_{0}(\beta)$
$\displaystyle=\prod_{i}\frac{1}{2\sinh\omega_{i}\beta/2}$
$\displaystyle=\prod_{i}e^{-\omega_{i}\tau/2}e^{-i\omega_{i}t/2}\frac{1}{1-e^{-\omega_{i}\beta}}$
$\displaystyle=Z^{\prime}_{0}(\beta)\det[e^{-\mathbf{W}\tau/2}e^{-i\mathbf{W}t/2}],$

where $Z_{0}^{\prime}$ is the partition function with the leading exponential
factor removed, i.e. with respect to the zero-point energy. The purely
quadratic contribution is then

$\displaystyle(*)$ $\displaystyle=\frac{(-1)^{n}}{Z^{\prime 2}_{0}}\det[\xi
e^{i\mathbf{W}t}]\det\left[\mathcal{L}\underline{e}\mathcal{L}^{-1}-\underline{\xi}^{-1}\right]^{-1}$
$\displaystyle=\frac{(-1)^{n}}{Z^{\prime
2}_{0}}e^{it\text{Tr}\mathbf{W}}\det\left[\left[\begin{array}[]{cc}\xi^{-1/2}&0\\\
0&1\end{array}\right]\mathcal{L}\left[\begin{array}[]{cc}e^{-}&0\\\
0&e^{+}\end{array}\right]\mathcal{L}^{-1}\left[\begin{array}[]{cc}\xi^{-1/2}&0\\\
0&1\end{array}\right]-\left[\begin{array}[]{cc}1&0\\\
0&\xi^{-1}\end{array}\right]\right]^{-1}$
$\displaystyle=\frac{(-1)^{n}}{Z^{\prime
2}_{0}}e^{it\text{Tr}\mathbf{W}}\det\left[\begin{array}[]{cc}\frac{1}{4}\xi^{-1/2}(\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{-}^{T})\xi^{-1/2}-1&-\text{transpose}\\\
\frac{1}{4}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{-}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{+}^{T})\xi^{-1/2}&\frac{1}{4}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T})-\xi^{-1}\end{array}\right]^{-1}$
$\displaystyle=\frac{1}{Z^{\prime
2}_{0}}e^{it\text{Tr}\mathbf{W}}\det\left[\begin{array}[]{cc}1-\frac{1}{4}\xi^{-1/2}(\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{-}^{T})\xi^{-1/2}&+\text{transpose}\\\
\frac{1}{4}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{-}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{+}^{T})\xi^{-1/2}&\frac{1}{4}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T})-\xi^{-1}\end{array}\right]^{-1}$

Before continuing, we can already see that in the low $T$ limit,
$\xi^{-1}\rightarrow 0$ and the determinant reduces to the vacuum expectation
value result.

We now factor this matrix determinant into that of the lower-right block and
its Schur complement,

$\displaystyle\det\left[\begin{array}[]{cc}A&B\\\ C&D\end{array}\right]$
$\displaystyle=\det(D)\det(A-BD^{-1}C).$

The first factor is

$\displaystyle\det\left[\frac{1}{4}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T})-\xi^{-1}\right]$
$\displaystyle=\det\left[\frac{1}{4}\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}(1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1}))\right]$
$\displaystyle=\det(\mathbf{\Lambda}_{+}/2)^{2}e^{it\text{Tr}[\mathbf{\Sigma\Omega}]}\det\left[1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})\right]$

As with the zero-temperature case, the last factor is calculated by
multiplying its eigenvalues, each of which does not cross the negative real
axis. The principal square root can then be taken without crossing branch-
cuts.

The second factor, the Schur complement, is

$\displaystyle\det\left[1-\frac{1}{4}\xi^{-1/2}\left(\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{-}^{T}\right.\right.$
$\displaystyle\qquad\left.\left.+(\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{-}^{T})(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{-}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{+}^{T})\right)\xi^{-1/2}\right].$

In the $e^{+}\rightarrow\infty$ limit, all diverging terms cancel exactly. To
remove all instances of $e^{+}$, we evaluate the total expression as

$\displaystyle\det[1-\frac{1}{4}\xi^{-1/2}(\cdots)\xi^{-1/2}],$

where

$\displaystyle(\cdots)$ $\displaystyle=(1)+(2)+(3),$ $\displaystyle(1)$
$\displaystyle=\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{+}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{-}^{T}-\mathbf{\Lambda}_{-}e^{+}\mathbf{\Lambda}_{-}^{T}$
$\displaystyle=\mathbf{\Lambda}_{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1})(1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1}))^{-1}\mathbf{\Lambda}^{\prime}_{+}\mathbf{\Lambda}_{-}^{T}$
$\displaystyle(2)$
$\displaystyle=\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{-}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{+}^{T}+\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{+}^{T}$
$\displaystyle=\mathbf{\Lambda}_{+}e^{-}\left[\mathbf{\Lambda}_{-}^{T}(1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1}))^{-1}\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}\mathbf{\Lambda}_{-}e^{-}+1\right]\mathbf{\Lambda}_{+}^{T}$
$\displaystyle(3)$
$\displaystyle=-\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{-}^{T}(\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{+}^{T}-\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}-4\xi^{-1})^{-1}\mathbf{\Lambda}_{+}e^{+}\mathbf{\Lambda}_{-}^{T}+\text{transpose}$
$\displaystyle=-\mathbf{\Lambda}_{+}e^{-}\mathbf{\Lambda}_{-}^{T}(1-\mathbf{\Lambda}^{\prime}_{+}e^{-}\mathbf{\Lambda}_{+}^{-1}(\mathbf{\Lambda}_{-}e^{-}\mathbf{\Lambda}_{-}^{T}+4\xi^{-1}))^{-1}\mathbf{\Lambda}^{\prime}_{+}\mathbf{\Lambda}_{-}^{T}+\text{transpose}$

We again assert that the eigenvalues of this matrix do not cross the negative
real axis, which provides for a branch-cut free factorization of the principal
square root.

## References

  * Sharp and Rosenstock [1964] T. E. Sharp and H. M. Rosenstock, “Franck-Condon Factors for Polyatomic Molecules,” J. Chem. Phys. 41, 3453 (1964).
  * Doktorov, Malkin, and Man’ko [1975] E. V. Doktorov, I. A. Malkin, and V. I. Man’ko, “Dynamical symmetry of vibronic transitions in polyatomic molecules and the Franck-Condon principle,” J. Mol. Spectrosc. 56, 1–20 (1975).
  * Caldwell and Gordon [1980] J. W. Caldwell and M. S. Gordon, “An approach to polyatomic Franck-Condon integrals: Application to the photoelectron spectrum of water,” J. Mol. Spectrosc. 84, 503–519 (1980).
  * Fernández and Tipping [1989] F. M. Fernández and R. H. Tipping, “Multidimensional harmonic oscillator matrix elements,” J. Chem. Phys. 91, 5505 (1989).
  * Dierksen [2005] M. Dierksen, “FCFAST, Ver. 1.0, Universität Münster, Münster, Germany,” (2005).
  * Rabidoux, Eijkhout, and Stanton [2016] S. M. Rabidoux, V. Eijkhout, and J. F. Stanton, “A Highly-Efficient Implementation of the Doktorov Recurrence Equations for Franck–Condon Calculations,” J. Chem. Theory Comput. 12, 728–739 (2016).
  * [7] S. L. Li and D. G. Truhlar, “FCBand 2017; http://comp.chem.umn.edu/fcband/,” .
  * Gozem and Krylov [2021] S. Gozem and A. I. Krylov, “The ezSpectra suite: An easy-to-use toolkit for spectroscopy modeling,” WIREs Comput. Mol. Sci. , e1546 (2021).
  * Dierksen and Grimme [2005] M. Dierksen and S. Grimme, “An efficient approach for the calculation of Franck-Condon integrals of large molecules,” J. Chem. Phys. 122, 244101 (2005).
  * Santoro  _et al._ [2007] F. Santoro, R. Improta, A. Lami, J. Bloino, and V. Barone, “Effective method to compute Franck-Condon integrals for optical spectra of large molecules in solution,” J. Chem. Phys. 126, 084509 (2007).
  * Beck  _et al._ [2000] M. Beck, A. Jäckle, G. Worth, and H.-D. Meyer, “The multiconfiguration time-dependent Hartree (MCTDH) method: a highly efficient algorithm for propagating wavepackets,” Phys. Rep. 324, 1–105 (2000).
  * Chen and Guo [1999] R. Chen and H. Guo, “The Chebyshev propagator for quantum systems,” Comp. Phys. Comm. 119, 19–31 (1999).
  * Vibók and Balint-Kurti [1992] Á. Vibók and G. G. Balint-Kurti, “Parametrization of complex absorbing potentials for time-dependent quantum dynamics,” J. Phys. Chem. 96, 8712–8719 (1992).
  * Leforestier and Wyatt [1998] C. Leforestier and R. E. Wyatt, “Optical potential for laser induced dissociation,” J. Chem. Phys. 78, 2334 (1998).
  * Muga  _et al._ [2004] J. Muga, J. Palao, B. Navarro, and I. Egusquiza, “Complex absorbing potentials,” Phys. Rep. 395, 357–426 (2004).
  * Luis, Bishop, and Kirtman [2004] J. M. Luis, D. M. Bishop, and B. Kirtman, “A different approach for calculating Franck-Condon factors including anharmonicity,” J. Chem. Phys. 120, 813 (2004).
  * Tannor and Heller [1982] D. J. Tannor and E. J. Heller, “Polyatomic Raman scattering for general harmonic potentials,” J. Chem. Phys. 77, 202 (1982).
  * Ianconescu and Pollak [2004] R. Ianconescu and E. Pollak, “Photoinduced Cooling of Polyatomic Molecules in an Electronically Excited State in the Presence of Dushinskii Rotations,” J. Phys. Chem. A 108, 7778–7784 (2004).
  * Baiardi, Bloino, and Barone [2013] A. Baiardi, J. Bloino, and V. Barone, “General Time Dependent Approach to Vibronic Spectroscopy Including Franck-Condon, Herzberg-Teller, and Duschinsky Effects,” J. Chem. Theory Comput. 9, 4097–4115 (2013).
  * Banerjee, Kröner, and Saalfrank [2012] S. Banerjee, D. Kröner, and P. Saalfrank, “Resonance Raman and vibronic absorption spectra with Duschinsky rotation from a time-dependent perspective: Application to $\beta$-carotene,” J. Chem. Phys. 137, 22A534 (2012).
  * Gelabert  _et al._ [2000] R. Gelabert, X. Giménez, M. Thoss, H. Wang, and W. H. Miller, “A Log-Derivative Formulation of the Prefactor for the Semiclassical Herman-Kluk Propagator,” J. Phys. Chem. A 104, 10321–10327 (2000).
  * Horváthy [2011] P. A. Horváthy, “The Maslov correction in the semiclassical Feynman integral,” Cent. Eur. J. Phys. 9, 1–12 (2011).
  * Gherib, Genin, and Ryabinkin [2022] R. Gherib, S. N. Genin, and I. G. Ryabinkin, “On the Importance of Well-Defined Thermal Correlation Functions in Simulating Vibronic Spectra,” J. Phys. Chem. A 126, 3947–3956 (2022).
  * Note [1] The procedure is simply to take small enough time-steps to reliably detect a numerical discontinuity in the correlation function.
  * Gilmore [2008] R. Gilmore,  _Lie Groups, Physics, and Geometry_ (Cambridge University Press, 2008).
  * Gilmore [1974] R. Gilmore, “Baker‐Campbell‐Hausdorff formulas,” J. Math. Phys. 15, 2090 (1974).
  * Note [2] 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.
  * Iso, Ohta, and Suyama [2018] S. Iso, H. Ohta, and T. Suyama, “Secular terms in Dyson series to all orders of perturbation,” Prog. Theor. Exp. Phys. 2018, 83A01 (2018).
  * Changala [2022] P. B. Changala, “NITROGEN v2.1.2, github.com/bchangala/nitrogen,” (2022).
  * Noga and Bartlett [1987] J. Noga and R. J. Bartlett, “The full CCSDT model for molecular electronic structure,” J. Chem. Phys. 86, 7041–7050 (1987).
  * Kendall, Dunning, and Harrison [1992] R. A. Kendall, T. H. Dunning, and R. J. Harrison, “Electron affinities of the first‐row atoms revisited. Systematic basis sets and wave functions,” J. Chem. Phys. 96, 6796 (1992).
  * Kucharski  _et al._ [2001] S. A. Kucharski, M. Włoch, M. Musiał, and R. J. Bartlett, “Coupled-cluster theory for excited electronic states: The full equation-of-motion coupled-cluster single, double, and triple excitation method,” J. Chem. Phys. 115, 8263–8266 (2001).
  * Matthews  _et al._ [2020] D. A. Matthews, L. Cheng, M. E. Harding, F. Lipparini, S. Stopkowicz, T.-C. Jagau, P. G. Szalay, J. Gauss, and J. F. Stanton, “Coupled-cluster techniques for computational chemistry: The CFOUR program package,” J. Chem. Phys. 152, 214108 (2020).
  * Holland and Lyman [1987] R. Holland and J. L. Lyman, “Ultraviolet absorption spectrum of molecular fluorine,” J. Quant. Spectrosc. Radiat. Transfer 38, 79–80 (1987).
  * Burkholder  _et al._ [2019] J. B. Burkholder, S. P. Sander, J. Abbatt, J. R. Barker, C. Cappa, J. D. Crounse, T. S. Dibble, R. E. Huie, C. E. Kolb, M. J. Kurylo, V. L. Orkin, C. J. Percival, D. M. Wilmouth, and P. H. Wine, “Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies, Evaluation No. 19, JPL Publication 19-5, Jet Propulsion Laboratory, Pasadena, CA, http://jpldataeval.jpl.nasa.gov,”  (2019).
  * Burkholder [1993] J. B. Burkholder, “Ultraviolet absorption spectrum of HOCl,” J. Geophys. Res. 98, 2963–2974 (1993).
  * Barnes, Sinha, and Michelsen [1998] R. J. Barnes, A. Sinha, and H. A. Michelsen, “Assessing the Contribution of the Lowest Triplet State to the Near-UV Absorption Spectrum of HOCl,” J. Phys. Chem. A 102, 8855–8859 (1998).
  * Minaev and Ågren [1998] B. F. Minaev and H. Ågren, “Response theory calculations of the singlet–triplet transition probabilities in the HOCl molecule,” J. Chem. Soc. Faraday Trans. 94, 2061–2067 (1998).
  * Pearson and Lovas [1977] R. Pearson and F. J. Lovas, “Microwave spectrum and molecular structure of methylenimine (CH2NH),” J. Chem. Phys. 66, 4149 (1977).
  * Bonačić-Koutecký and Michl [1985] V. Bonačić-Koutecký and J. Michl, “Photochemical syn-anti isomerization of a Schiff base: A two-dimensional description of a conical intersection in formaldimine,” Theor. Chim. Acta 68, 45–55 (1985).
  * Hazra, Chang, and Nooijen [2004] A. Hazra, H. H. Chang, and M. Nooijen, “First principles simulation of the UV absorption spectrum of ethylene using the vertical Franck-Condon approach,” J. Chem. Phys. 121, 2125 (2004).
  * Teslja, Nizamov, and Dagdigian [2004] A. Teslja, B. Nizamov, and P. J. Dagdigian, “The Electronic Spectrum of Methyleneimine,” J. Phys. Chem. A 108, 4433–4439 (2004).
  * Matthews and Stanton [2016] D. A. Matthews and J. F. Stanton, “A new approach to approximate equation-of-motion coupled cluster with triple excitations,” J. Chem. Phys. 145, 124102 (2016).
  * Ndengué  _et al._ [2021] S. Ndengué, E. Quintas-Sánchez, R. Dawes, and D. Osborn, “The Low-Lying Electronic States of NO2: Potential Energy and Dipole Surfaces, Bound States, and Electronic Absorption Spectrum,” J. Phys. Chem. A 125, 5519–5533 (2021).
  * Bogumil  _et al._ [2003] K. Bogumil, J. Orphal, T. Homann, S. Voigt, P. Spietz, O. C. Fleischmann, A. Vogel, M. Hartmann, H. Kromminga, H. Bovensmann, J. Frerick, and J. P. Burrows, “Measurements of molecular absorption spectra with the SCIAMACHY pre-flight model: instrument characterization and reference data for atmospheric remote-sensing in the 230–2380 nm region,” J. Photochem. Photobiol. A 157, 167–184 (2003).
  * Keller-Rudek  _et al._ [2013] H. Keller-Rudek, G. K. Moortgat, R. Sander, and R. Sörensen, “The MPI-Mainz UV/VIS spectral atlas of gaseous molecules of atmospheric interest,” Earth Syst. Sci. Data 5, 365–373 (2013).
  * Note [3] In fact, the use of a vertical force field, rather than that of the upper-state minimum, is in some sense a crude form of renormalization.
  * Ruscic [1986] B. Ruscic, “Fourier transform photoelectron spectroscopy: The correlation function and the harmonic oscillator approximation,” J. Chem. Phys. 85, 3776–3784 (1986), https://doi.org/10.1063/1.450898 .
  * Siebrand [1967] W. Siebrand, “Radiationless Transitions in Polyatomic Molecules. I. Calculation of Franck—Condon Factors,” J. Chem. Phys. 46, 440 (1967).
  * Fischer [2003] S. Fischer, “Correlation Function Approach to Radiationless Transitions,” J. Chem. Phys. 53, 3195 (2003).
  * Peng  _et al._ [2007] Q. Peng, Y. Yi, Z. Shuai, and J. Shao, “Excited state radiationless decay process with Duschinsky rotation effect: Formalism and implementation,” J. Chem. Phys. 126, 114302 (2007).
  * Hougen and Watson [1965] J. T. Hougen and J. K. G. Watson, “Anomalous rotational line intensities in electronic transitions of polyatomic molecules: axis-switching,” Can. J. Phys. 43, 298–320 (1965).