On-the-fly ab initio three thawed Gaussians approximation: A semiclassical approach to Herzberg-Teller spectra

TL;DR

This paper introduces a novel on-the-fly ab initio semiclassical method using three thawed Gaussians to improve the accuracy of vibronic spectra calculations, especially for weakly-allowed transitions.

Contribution

The paper develops and applies a three thawed Gaussians approximation to better capture anharmonic effects in vibronic spectra beyond traditional harmonic models.

Findings

Improved spectral accuracy over harmonic and Condon models.

Effective inclusion of Herzberg-Teller effects.

Enhanced modeling of vibrationally resolved spectra.

Abstract

Evaluation of symmetry-forbidden or weakly-allowed vibronic spectra requires treating the transition dipole moment beyond the Condon approximation. To include vibronic spectral features not captured in the global harmonic models, we have recently implemented an on-the-fly ab initio extended thawed Gaussian approximation, where the propagated wavepacket is a Gaussian multiplied by a linear polynomial. To include more anharmonic effects, here we represent the initial wavepacket by a superposition of three independent Gaussian wavepackets - one for the Condon term and two displaced Gaussians for the Herzberg-Teller part. Application of this ab initio "three thawed Gaussians approximation" to vibrationally resolved electronic spectra of the phenyl radical and benzene shows a clear improvement over the global harmonic and Condon approximations. The orientational averaging of spectra, the relation between the gradient of the transition dipole moment and nonadiabatic coupling vectors, and the details of the extended and three thawed Gaussians approximation are discussed.

Figures

Equations

$$\sigma (\overrightarrow{ϵ},\omega )=\frac{2\pi \omega }{\hslash c}{\int }_{-\infty }^{\infty }{C}_{\mu \mu }(\overrightarrow{ϵ},t)e^{i\omega t}dt$$

$$C_{\mu \mu }(\overrightarrow{ϵ},t)=⟨1,g|e^{i{\hat{H}}_{1}t/\hslash }{\hat{\mu }}_{12}{e}^{-i{\hat{H}}_{2}t/\hslash }{\hat{\mu }}_{21}|1,g⟩,$$

$${\hat{\mu }}_{21}:={\hat{\vec{\mu }}}_{21}\cdot \overrightarrow{ϵ}$$

$$\sigma (\overrightarrow{ϵ},\omega )={\overrightarrow{ϵ}}^T\cdot \overleftrightarrow{\sigma }(\omega )\cdot \overrightarrow{ϵ}.$$

$$\overline{\sigma (\overrightarrow{ϵ},\omega )}$$

$$=\frac{1}{3}\left[{\sigma }_{xx}(\omega )+{\sigma }_{yy}(\omega )+{\sigma }_{zz}(\omega )\right]$$

$$=\frac{1}{3}\left[\sigma ({\vec{e}}_x,\omega )+\sigma ({\vec{e}}_y,\omega )+\sigma ({\vec{e}}_z,\omega )\right].$$

$$\overline{C_{\mu \mu }(\overrightarrow{ϵ},t)}=\frac{1}{3}\mathrm{Tr}[{\overleftrightarrow{C}}_{\mu \mu }(t)].$$

$$\overline{C_{\mu \mu }(\overrightarrow{ϵ},t)}=\frac{1}{3}\left[C_{{\mu }_x{\mu }_x}(t)+C_{{\mu }_y{\mu }_y}(t)+C_{{\mu }_z{\mu }_z}(t)\right],$$

$$\overline{C_{\mu \mu }(\overrightarrow{ϵ},t)}=\frac{1}{3}{C}_{|\vec{\mu }||\vec{\mu }|}(t),$$

$$\vec{\mathbf{\mu }}(q)\approx \vec{\mathbf{\mu }}(q_{\text{eq}}).$$

$$\vec{\mathbf{\mu }}(q)\approx \vec{\mathbf{\mu }}(q_{\text{eq}})+{\partial }_q\vec{\mathbf{\mu }}{|}_{q_{\text{eq}}}^T\cdot (q-q_{\text{eq}}),$$

$${\partial }_j\vec{\mathbf{\mu }}(q)=\left[\vec{\mathbf{\mu }}(q),{\mathbf{F}}_j(q)\right]+{\partial }_j{\vec{\mu }}_{\text{nu}}(q)1$$

$${\vec{\mu }}_{\alpha \beta }(q):=⟨\alpha (q)|{\hat{\vec{\mu }}}_{\text{mol}}|\beta (q)⟩,$$

$$F_{\alpha \beta ,j}:=⟨\alpha (q)|{\partial }_j\beta (q)⟩,$$

$${\vec{\mu }}_{\text{nu}}(q)=e\sum _{N=1}^{N_{\text{nu}}}{Z}_N{\vec{R}}_N(q).$$

$${\partial }_j{\vec{\mu }}_{\alpha 1}\approx -\sum _{\gamma }{F}_{\alpha \gamma ,j}{\vec{\mu }}_{\gamma 1}$$

$$C_{\mu \mu }(\overrightarrow{ϵ},t)=C(t)e^{i{E}_{1,g}t/\hslash }$$

$$C(t)=⟨\varphi (0)|\varphi (t)⟩$$

$$|\varphi (t)⟩=e^{-i{\hat{H}}_{2}t/\hslash }|\varphi (0)⟩.$$

$${\psi }_t(q)=N_0\exp \left\{-(q-q_t{)}^T\cdot A_t\cdot (q-q_t)+\frac{i}{\hslash }\left[p_t^T\cdot (q-q_t)+{\gamma }_t\right]\right\}$$

$${\hat{H}}_{\text{eff}}(t)\equiv H_{\text{eff}}(\hat{q},\hat{p},t)=\frac{1}{2}{\hat{p}}^T\cdot m^{-1}\cdot \hat{p}+V_{\text{eff}}(\hat{q},t),$$

$$V_{\text{eff}}(q,t)=V{|}_{q_t}+({\partial }_{q}V{|}_{q_t}{)}^T\cdot (q-q_t)+\frac{1}{2}(q-q_t{)}^T\cdot {\mathrm{Hess}}_{q}V{|}_{q_t}\cdot (q-q_t).$$

$${\left({\mathrm{Hess}}_{q}V\right)}_{ij}:={\partial }_i{\partial }_{j}V.$$

$${\dot{q}}_t$$

$${\dot{p}}_t$$

$${\dot{A}}_t$$

$${\dot{\gamma }}_t$$

$$M_t=\left(\begin{array}{cc}{M}_{t,qq}& M_{t,qp}\\ M_{t,pq}& M_{t,pp}\end{array}\right):=\left(\begin{array}{cc}\frac{\partial q_t}{\partial q_0}& \frac{\partial q_t}{\partial p_0}\\ \frac{\partial p_t}{\partial q_0}& \frac{\partial p_t}{\partial p_0}\end{array}\right),$$

$${\varphi }_0(q)=P\left(q-q_0\right){\psi }_0(q).$$

$${\varphi }_0(q)=P\left(\frac{\hslash }{i}\frac{\partial }{\partial p_0}\right){\psi }_0(q).$$

$${\varphi }_t(q)=P\left(\frac{\hslash }{i}\frac{\partial }{\partial p_0}\right){\psi }_t(q),$$

$${\varphi }_t(q)=\left[{\mu }_{21}(q_0)+b_t^T\cdot (q-q_t)\right]{\psi }_t(q),$$

Full text

I Introduction

Electronic spectroscopy lies at the core of modern physical chemistry; not only has it been the driving force in developing first insights into the atomic and molecular structure Quack and Merkt (2011); Herzberg (1966), but it is also the method of choice for unraveling essential chemical and physical processes. The time-dependent approach to electronic spectroscopy obtains the spectrum as a Fourier transform of an appropriate time correlation function Heller (1981a); Mukamel (1995), which requires performing molecular quantum dynamics, but provides more dynamical information about the interaction of the molecule with light. Even though only relatively short (femtosecond to picosecond) time scales are typically involved in electronic spectroscopy, exact methods for solving the time-dependent Schrödinger equation become quickly impractical for larger molecules due to the exponential scaling with the number of degrees of freedom Meyer, Gatti, and Worth (2009); Gatti (2014). In contrast, approximate semiclassical methods reduce the exponential quantum problem to the propagation of classical trajectories, while often maintaining sufficient accuracy, at least at short times. With the advent of computationally feasible, though not necessarily cheap, ab initio calculations, semiclassical dynamics benefits from yet another advantage over quantum dynamical methods: the formulation in terms of classical trajectories allows for an efficient on-the-fly implementation, overcoming the need for exploring the full potential energy surface.

Among many semiclassical approximations Miller (1970, 2001); Herman and Kluk (1984); Kay (2005); Zhang and Pollak (2004); Grossmann (2006); Mollica and Vaníček (2011); Ceotto, Di Liberto, and Conte (2017); Buchholz, Grossmann, and Ceotto (2018), one of the simplest, yet very intuitive and often surprisingly accurate, is the thawed Gaussian approximation (TGA) of Heller Heller (1975). Although his work on Gaussian wavepackets Heller (1981b) inspired a number of “direct”, i.e., on-the-fly ab initio quantum Ben-Nun, Quenneville, and Martínez (2000); Saita and Shalashilin (2012); Makhov et al. (2017); Richings et al. (2015); Curchod and Martínez (2018) and semiclassical Tatchen and Pollak (2009); Ceotto et al. (2009a, b); Wong et al. (2011); Ianconescu, Tatchen, and Pollak (2013); Gabas, Conte, and Ceotto (2017) methods, the original thawed Gaussian approximation was only revisited recently, when it was implemented within the on-the-fly ab initio framework for evaluating molecular absorption, emission, and photoelectron spectra Wehrle, Šulc, and Vaníček (2014); Wehrle, Oberli, and Vaníček (2015). Due to its computational efficiency in comparison with exact quantum dynamics methods, the TGA can account for Wehrle, Šulc, and Vaníček (2014); Wehrle, Oberli, and Vaníček (2015) the full dimensionality of the system and mode-mixing (Duschinsky effect), both of which can significantly influence the wavepacket dynamics Mahapatra et al. (2005). Recently, we have implemented Patoz, Begušić, and Vaníček (2018) an extension to the on-the-fly ab initio thawed Gaussian approximation, aimed at describing the Herzberg–Teller contribution to the vibrational structure of electronic absorption spectra. The extended thawed Gaussian approximation (ETGA) Lee and Heller (1982); Patoz, Begušić, and Vaníček (2018), propagates a Gaussian wavepacket multiplied by a general polynomial in nuclear coordinates using a local harmonic approximation, i.e., the second-order expansion of the potential about the center of the wavepacket. Such a wavepacket, however, is wider than a simple Gaussian wavepacket, which raises the question whether more than a single guiding classical trajectory should be used to correctly account for the anharmonicity of the potential energy surface.

Here, we present the on-the-fly ab initio implementation of an alternative method for treating Herzberg–Teller spectra beyond the Condon approximation. This semiclassical method is also based on the thawed Gaussian approximation, yet, unlike the extended TGA, resolves the initial Herzberg–Teller wavepacket into three well-defined Gaussians and propagates them independently. This “three thawed Gaussians approximation” (3TGA) is applied to evaluate the absorption spectra of the phenyl radical and benzene. The 3TGA is compared with the Condon approximation in order to analyze the importance of the Herzberg–Teller contribution, with the global harmonic approaches to assess the extent of anharmonicity, and with the extended TGA to analyze the effects of wavepacket splitting.

Notation: As the analysis of electronic spectra involves three different vector spaces, let us summarize the notation for reference here, although most should be clear from the context. If $D$ denotes the number of nuclear degrees of freedom and $S$ the number of electronic states considered, the three vector spaces are the nuclear $D$-dimensional real coordinate space $\mathbb{R}^{D}$, electronic $S$-dimensional complex Hilbert space $\mathbb{C}^{S}$, and the ambient $3$-dimensional space $\mathbb{R}^{3}$. To distinguish these spaces, $3$-dimensional vectors will be denoted with an arrow (e.g., $\vec{\epsilon}$), whereas $D$-dimensional nuclear vectors or $D\times D$ matrices will use no special notation (e.g., generalized nuclear coordinates $q$). Scalar and matrix products in both the $3$-dimensional and $D$-dimensional nuclear spaces will be denoted with a dot (as in $\vec{\mu}_{21}\cdot\vec{\epsilon}$ or $p^{T}\cdot m^{-1}\cdot p$). We shall use the bold font for $S\times S$ matrices (such as $\bm{\mu}$) representing electronic operators expressed in the $S$-state basis of the electronic Hilbert space. The matrix product will use no special notation; it will be expressed by a juxtaposition of the matrices (as in $\mathbf{AB}$).

II Theory

II.1 Time-dependent approach to spectroscopy

Within the electric-dipole approximation and first-order time-dependent perturbation theory, the absorption cross section for a linearly polarized light of frequency $\omega$ can be expressed as the Fourier transform

[TABLE]

of the dipole time autocorrelation function $C_{\mu\mu}(\vec{\epsilon},t)$. Assuming the zero temperature approximation, the initial state is $|1,g\rangle$, i.e., the ground vibrational state $g$ of the ground electronic state $1$. Assuming, furthermore, that the incident radiation is in resonance only with a single pair of electronic states $1$ and $2$ that are not vibronically coupled, the dipole time autocorrelation function reduces to Heller (1981a)

[TABLE]

where $\hat{H}_{1}$ and $\hat{H}_{2}$ are nuclear Hamiltonian operators in the ground and excited electronic states, and

[TABLE]

is the projection of the matrix element $\hat{\vec{\mu}}_{21}$ of the molecular transition dipole moment matrix $\hat{\vec{\bm{\mu}}}$ on a three-dimensional polarization unit vector $\vec{\epsilon}$ of the electric field.

II.2 Orientational averaging of the

spectrum

In gas phase or another isotropic medium, one has to average the vibronic spectrum (1) over all orientations of the molecule with respect to the polarization $\vec{\epsilon}$ of the electric field. Within the Condon approximation, the transition dipole moment is independent of coordinates, and this averaging is trivial, but for a general dipole moment, one has to be more careful. It is surprising that in theoretical papers on vibronic spectroscopy, the averaging is often ignored despite an extensive work on orientational averaging of both linear and nonlinear spectra Andrews and Thirunamachandran (1977); Gelin, Borrelli, and Domcke (2017).

It is useful to define the spectrum tensor $\overleftrightarrow{\sigma}(\omega)$, from which the spectrum (1) for a specific polarization $\vec{\epsilon}$, is obtained by “evaluation”:

[TABLE]

The orientational averaging of the spectrum corresponds to the averaging of $\sigma(\vec{\epsilon},\omega)$ over all unit vectors $\vec{\epsilon}$. Remarkably, due to the isotropy of the 3-dimensional Euclidean space, the average over all orientations does not have to be performed numerically but can be reduced to a simple arithmetic average over only three arbitrary orthogonal orientations of the molecule with respect to the field:

[TABLE]

Appendix A contains an explicit proof of these equalities, of which the first is coordinate-independent and the other two are explicit in Cartesian coordinates ($\vec{e}_{x}$ denotes the unit vector along the $x$-axis.). This final result for the orientational average is exact for any coordinate dependence of the transition dipole operator and also for any quantum or semiclassical dynamical method for evaluating the autocorrelation function. Moreover, it is valid even for arbitrary pure or mixed initial molecular states and arbitrary nonadiabatic or electric-dipole couplings among the $S$ electronic states. As the Fourier transform is a linear operation, one may choose to apply the averaging already to the dipole autocorrelation function instead of the spectrum:

[TABLE]

Now let us go back to our two-state system, where only one element of the electric transition dipole moment, $\vec{\mu}_{21}(q)$ plays a role and $C_{\mu\mu}(\vec{\epsilon},t)$ is given by Eq. (2). Then it is convenient to suppress the subscripts $21$ and denote the $21$ element simply as $\vec{\mu}(q)$, which we shall do for the remainder of this subsection. Among a continuum of other possibilities, expression (6) provides two simple yet exact recipes for the orientational average: one can take the arithmetic average of $C_{\mu\mu}(\vec{\epsilon},t)$ either for three orthogonal polarizations $\vec{\epsilon}$ and a fixed molecular orientation Hein et al. (2012),

[TABLE]

or for three orthogonal orientations of the molecule and a fixed polarization $\vec{\epsilon}$. Within the Condon approximation, in which $\vec{\mu}$ is coordinate-independent, the orientational average (7) simplifies further into the standard textbook recipe

[TABLE]

where $|\vec{\mu}|$ is the magnitude of the transition dipole moment, so only a single calculation is required—for the transition dipole moment aligned with the field.

II.3 Condon and Herzberg–Teller approximations for

the transition dipole moment

The electric transition dipole moment is, in general, a function of nuclear coordinates, yet, within the Condon approximation Condon (1926), this moment is assumed to be independent of the molecular geometry and is approximated to the zeroth order around the initial geometry:

[TABLE]

The widespread use of this approximation is justified by its validity in many systems, as it can describe most of the strongly symmetry-allowed transitions both qualitatively and quantitatively. However, a number of molecules exhibit symmetry-forbidden (also called “electronically forbidden”) transitions, i.e., transitions $\alpha\leftarrow\beta$ with $\vec{\mu}_{\alpha\beta}(q_{\text{eq}})=0$, which cannot be described within the Condon approximation. Such systems, as well as systems in which the Condon term is small but not exactly zero, can be treated with the Herzberg–Teller approximation Herzberg and Teller (1933) that takes into account also the gradient of the transition dipole moment with respect to nuclear degrees of freedom:

[TABLE]

where $\partial_{q}$ is an abbreviation for $\partial/\partial q$.

II.4 Relationship between the gradient of the

transition dipole moment and nonadiabatic vector couplings

Although the Herzberg–Teller approximation is often discussed in terms of vibronic couplings between different electronic states Seidner et al. (1992), this relation is not obvious from Eq. (10). Here we demonstrate this relationship explicitly. In particular, we shall prove a remarkable equality

[TABLE]

satisfied by the gradient of the matrix representation $\vec{\bm{\mu}}(q)$ of the molecular dipole operator $\hat{\vec{\mu}}_{\text{mol}}$ at the nuclear configuration $q$. In this equation, matrix elements of $\vec{\bm{\mu}}(q)$ are defined as partial, electronic expectation values

[TABLE]

elements of the matrix $\mathbf{F}_{j}(q)$ of nonadiabatic vector couplings are obtained as

[TABLE]

and $\vec{\mu}_{\text{nu}}$ is the nuclear component of $\hat{\vec{\mu}}_{\text{mol}}$:

[TABLE]

Here $Z_{N}$ is the atomic number and $\vec{R}_{N}(q)$ the Cartesian coordinates of the $N$th nucleus expressed as a function of the generalized nuclear coordinates $q$, which can be Cartesian, normal-mode, Jacobi, or any other coordinates.

Direct interpretation of relation (11), which is proven in Appendix B, explains the concepts of vibronic transitions and intensity borrowing. Namely, the gradient of the transition dipole moment between states $\alpha$ and $\beta$ can be nonzero only if there exists an intermediate state $\gamma$ that is nonadiabatically (i.e., vibronically) coupled to one of the states and electric-dipole coupled to the other state. This is even more explicit in the component expression (50) in the penultimate line of the proof in the Appendix. Typically, the nonadiabatic couplings with the ground state ($\beta=1$) can be neglected at the ground-state optimized geometry, around which the transition dipole moment is expanded. This leads to an expression Li et al. (2010)

[TABLE]

that explains the meaning of intensity borrowing, in which the symmetry forbidden transition to the formally dark state $\alpha$ occurs due to a borrowing of the transition intensity from the neighboring bright electronic states $\gamma$ that are vibronically coupled to the state $\alpha$. Note that, despite introducing nonadiabatic couplings between excited electronic states, the original Born-Oppenheimer picture adopted in Subsection II.1 remains valid; the vibronic couplings that induce the transition do not necessarily influence the nuclear wavepacket dynamics. Although these nonadiabatic couplings are essential for describing the existence of symmetry-forbidden spectra, they are still rather small and their contribution to the field-free Hamiltonian of the system is negligible. The rather high resolution of the absorption spectrum of benzene supports these considerations—otherwise, significant population transfer would lead to the shortening of the excited-state lifetime and, consequently, to significant broadening of the spectral lines.

II.5 Extended thawed Gaussian approximation

To evaluate the dipole time autocorrelation function (2) corresponding to a specific polarization of the electric field, let us rewrite it as

[TABLE]

in terms of the vibrational zero-point energy $E_{1,g}$ of the ground electronic state and the wavepacket autocorrelation function

[TABLE]

for the initial wavepacket $|\phi(0)\rangle=\hat{\mu}_{21}|1,g\rangle$ propagated on the excited-state potential energy surface with the Hamiltonian $\hat{H}_{2}$:

[TABLE]

This propagation is the most demanding part of the spectra calculation and we shall use an on-the-fly ab initio semiclassical method for this purpose. One of the simplest semiclassical methods, the TGA Heller (1975), relies on the fact that a Gaussian wavepacket conserves its form when propagated in a time-dependent harmonic potential; it extends the global harmonic methods by approximating the potential energy surface locally to the second order, in what is known as the local harmonic approximation. Although approximating the initial state in electronic spectroscopy with a Gaussian is only reasonable within the Condon approximation, let us first discuss this simplest case because it will serve as a starting point for extensions to more general forms of the initial wavepacket needed to describe Herzberg–Teller spectra.

Gaussian wavepacket considered in TGA is given in the position representation as

[TABLE]

where $N_{0}$ is a normalization constant, $(q_{t},p_{t})$ are the phase-space coordinates of the center of the wavepacket, $A_{t}$ is a symmetric complex width matrix, and $\gamma_{t}$ a complex number whose real part is a dynamical phase and imaginary part guarantees the normalization for times $t>0$. The wavepacket is propagated with a Hamiltonian

[TABLE]

where $m$ is the diagonal mass matrix and $V_{\text{eff}}$ an effective time-dependent potential given by the local harmonic approximation of the true potential $V$ at the center of the wavepacket:

[TABLE]

To avoid the confusion between the Laplacian and Hessian operators, both of which are often denoted by $\partial^{2}/\partial q^{2}$, we use the notation $\operatorname{Hess}_{q}V$ for Hessian matrix with respect to nuclear coordinates:

[TABLE]

Insertion of the wavepacket ansatz (19) and the effective potential (21) into the time-dependent Schrödinger equation yields a system of ordinary differential equations Heller (1975)

[TABLE]

for evolving parameters of the Gaussian; $L_{t}$ denotes the Lagrangian. This system of differential equations, which is within the local harmonic approximation equivalent to the solution of the Schrödinger equation, implies that phase-space coordinates $q_{t}$ and $p_{t}$ follow classical equations of motion, while the propagation of the width $A_{t}$ and phase $\gamma_{t}$ involves propagating the $2D\times 2D$ stability matrix

[TABLE]

which depends on the Hessians of the potential energy surface $V$. See Appendix C for details.

To tackle a more general coordinate dependence of the transition dipole moment, such as the Herzberg–Teller approximation (10), Lee and Heller Lee and Heller (1982) proposed the extended TGA (ETGA) that considers a more general form of the initial wavepacket, namely a Gaussian multiplied by a polynomial Lee and Heller (1982); Patoz, Begušić, and Vaníček (2018):

[TABLE]

Otherwise, the ETGA is the same as the original TGA; it uses the local harmonic approximation (21) for the potential along the trajectory $q_{t}$, but makes no other approximation.

Because the only dependence of $\psi_{0}(q)$ on $p_{0}$ comes from the exponent $p_{0}^{T}\cdot\left(q-q_{0}\right)$ [see Eq. (19)], the polynomial prefactor in Eq. (27) can be replaced by the same polynomial in the derivatives with respect to $p_{0}$:

[TABLE]

In Appendix D, we prove that the local harmonic approximation implies that the ETGA wavepacket at time $t$ has the simple form Lee and Heller (1982)

[TABLE]

where the dependence of $\psi_{t}(q)$ on initial conditions $q_{0}$ and $p_{0}$ is taken into account. In particular, equations of motion (22)-(25) for $q_{t}$, $p_{t}$, $A_{t}$, and $\gamma_{t}$ remain unchanged.

As for the parameters of the polynomial, here we consider only the constant and linear terms required in the Herzberg–Teller approximation (10). In Appendix D, we prove that Lee and Heller (1982)

[TABLE]

where the linear parameter of the polynomial at time $t$ is

[TABLE]

Because all ingredients needed in Eq. (31) are already evaluated for the propagation of the parameters of the Gaussian, the evaluation of $b_{t}$ comes at almost no additional cost.

II.6 Three thawed Gaussians approximation

Let us now describe another generalization of the TGA, which, as the ETGA, can also evaluate the Herzberg–Teller spectra, but, in contrast to the ETGA, can also account for wavepacket splitting. Recall that the ETGA wavepacket is guided by a single trajectory propagated using the local harmonic approximation of the potential at the center of the wavepacket. Whereas the center of a Gaussian maximizes the probability density, this is false for purely Herzberg–Teller wavepackets [Eq. (30) with $\mu_{21}(q_{0})=0$], where the probability density at the center is zero. Although the ETGA is still an exact solution of the time-dependent Schrödinger equation in a global harmonic potential, its performance in a general potential can, therefore, be questioned Lee and Heller (1982). To further investigate the errors introduced by the local harmonic approximation, we propose to approximate the Herzberg–Teller part

[TABLE]

of the initial state used in the ETGA with an antisymmetric linear combination of two displaced Gaussians (see Fig. 1):

[TABLE]

where $g_{q_{c}}(q)$ is a normalized Gaussian $\psi_{0}(q)$ centered at $q_{c}$ instead of $q_{0}$, with zero initial momentum and phase ($p_{0}=\gamma_{0}=0$):

[TABLE]

$\Delta_{d}$ is a displacement vector and $f_{d}$ a scaling factor. Note also that the same width matrix $A_{0}$ is used for $g_{q_{c}}$ as for $\psi_{0}(q)$.

We will now show that the approximative wavepacket $\phi_{0}^{\text{3TGA-HT}}$ converges to the exact Herzberg–Teller wavepacket $\phi_{0}^{\text{ETGA-HT}}$when the displacement $\Delta_{d}$ approaches zero along an appropriate direction and a corresponding scaling factor $f_{d}$ is used. For brevity, throughout this subsection we shall again write $\mu$ instead of $\vec{\mu}_{21}(q_{0})\cdot\vec{\epsilon}$ and $\partial_{q}\mu$ instead of $\partial_{q}(\vec{\mu}_{21}(q)\cdot\vec{\epsilon})|_{q_{0}}$.

First, let us discuss the shape of the exact wavepacket $\phi_{0}^{\text{ETGA-HT}}$. Obviously, in one nuclear dimension, $\phi_{0}^{\text{ETGA-HT}}$ has exactly two local extrema: a minimum and a maximum. In Appendix E we prove that this property holds for any number $D$ of nuclear degrees of freedom: there are two and only two local extrema, a maximum and minimum, located at

[TABLE]

with the displacement vector

[TABLE]

It is rather obvious that in order for $\phi_{0}^{\text{3TGA-HT}}$ to be a good fit of $\phi_{0}^{\text{ETGA-HT}}$, the two Gaussians in Eq. (33) must be displaced from $q_{0}$ along the direction of the extrema of $\phi_{0}^{\text{ETGA-HT}}$. We therefore choose

[TABLE]

where the dimensionless parameter $d$ controls the magnitude of the displacement. The scaling factor $f_{d}$ is obtained by equating the norms of $\phi_{0}^{\text{3TGA-HT}}$ and $\phi_{0}^{\text{ETGA-HT}}$, which gives

[TABLE]

A quick inspection reveals $\phi_{0}^{\text{ETGA-HT}}$ to be a directional derivative of $g_{q_{0}}(q)$ and $\phi_{0}^{\text{3TGA-HT}}$ a finite-difference approximation of this derivative, exact in the limit $d\rightarrow 0$. One can see this explicitly by noting that the normalized initial overlap of the two wavepackets,

[TABLE]

is maximized and converges to unity as $d$ approaches zero. Either proof implies that a better description of the initial wavepacket is obtained with smaller values of $d$. On the other hand, anharmonicity effects are included only when the two initial wavepackets are significantly displaced, implying that a larger $d$ is needed. A natural displacement that is neither too small nor too large corresponds to placing the two Gaussians at the extrema of the wavepacket by setting $d=1$ and $\Delta_{d}=\Delta q$; in this case, the initial overlap Eq. (39) is $\approx 98\%$. As this choice of $d$ results in a large initial overlap and in trajectories that correspond to the maxima of the probability density of the Herzberg–Teller component, we recommend using $d=1$ in all applications. However, in the calculations presented in Section IV, we will test not only $d=1$ but also several smaller displacements, giving even larger initial overlaps.

Finally, by adding the Condon term, one obtains a superposition of three Gaussian wavepackets that are propagated independently; the total initial wavepacket in the “three thawed Gaussians approximation” (3TGA) is

[TABLE]

In contrast to both the original and extended TGA, where the evolution given by the time-dependent effective Hamiltonian is exactly unitary, the norm of the 3TGA wavepacket is not conserved; since the three Gaussians are propagated with three different effective Hamiltonians, their overlaps are time-dependent and change the norm of the total 3TGA wavepacket. To see this, consider several wavepackets $|\psi_{i}\rangle$, each propagated with its own evolution operator $\hat{U}_{i}$. If $\hat{U}_{i}\neq\hat{U}_{j}$, then the overlap of two such wavepackets is time-dependent:

[TABLE]

The time dependence of the norm of the 3TGA wavepacket arises from these time-dependent overlaps since

[TABLE]

Due to the time dependence of the norm, the calculated dipole time autocorrelation function has to be renormalized at each time step.

II.7 Ab initio implementation

The ab initio implementation of the thawed Gaussian approximation has been discussed in Refs. 29; 30; only a brief overview is given here. The propagation of a thawed Gaussian wavepacket requires a single classical trajectory on the excited-state potential energy surface and the Hessians of the excited-state potential along this trajectory. These data are evaluated in Cartesian coordinates and are then transformed to mass-scaled normal mode coordinates, after removing the translational degrees of freedom by translation to the center of mass frame and rotational degrees of freedom by rotation to the Eckart frame. The choice of the coordinates and correct coordinate transformation is the essence of the on-the-fly ab initio implementation; mass-scaled normal modes of the ground electronic state are useful not only for a straightforward construction of the initial wavepacket, but also for further interpretation of the results.

Within the 3TGA, each thawed Gaussian is propagated using the above mentioned scheme. The initial positions are obtained in the following way: First, the gradient of the transition dipole moment in Cartesian coordinates is projected onto one of three orthogonal polarizations of the electric field and transformed to the ground-state normal mode coordinates. Then the centers of the displaced Gaussians are found and transformed back to Cartesian coordinates for an on-the-fly ab initio propagation. In addition to the central trajectory, each polarization of the electric field requires, in general, two different displaced trajectories to be evaluated.

III Computational details

All ab initio calculations were performed using B3LYP/6-31+G(d,p) density functional theory for the ground state and time-dependent density functional theory for the excited state, as implemented in the Gaussian09 package Frisch et al. . The gradient of the transition dipole moment was computed using a finite-difference approach, as reported in Ref. 32. The ground-state potential energy surface was approximated with a global harmonic potential in order to obtain the initial vibrational state.

The orientational averaging of the 3TGA spectra requires additional trajectories because the gradient of the transition dipole moment depends on the orientation of the molecule. In general, one has to compute the spectrum for three orthogonal orientations of the molecule with respect to the electric field, which implies one central trajectory for the Condon contribution and two additional trajectories for each orientation to represent the Herzberg–Teller contribution to the wavepacket; seven trajectories in total. In phenyl radical, all seven trajectories have to be evaluated, while in benzene, since the transition is symmetry-forbidden and since the gradient of the $z$ component of the transition dipole moment is zero, only four trajectories are needed (two for $x$ and two for $y$ polarization). In contrast, the orientational averaging of the ETGA spectra requires only a single ab initio trajectory because the averaging can be performed by changing only the polynomial part of the wavepacket, which depends on the polarization of the electric field, but does not need additional ab initio data.

A time step of 8 a.u. ($\approx 0.194$ fs) was chosen for all trajectories; 3000 steps (giving a total time of 580 fs) were run for the phenyl radical and 5000 steps (970 fs) for benzene. The Hessians of the potential were computed every four steps, while the intermediate Hessians were obtained by second-order interpolation.

Gaussian broadening of the resulting spectra was used for the phenyl radical (half-width at half-maximum of 100 cm*-1*), while for benzene the autocorrelation function is multiplied by $\cos^{2}(\pi t/2T)$ ($T$ is the length of the simulation) because this function preserves most of the autocorrelation function and introduces only slight broadening. However, longer propagation would be required to simulate the spectrum with peaks as narrow as in the experiment. For comparison with experiment, unless otherwise stated, we scale and shift the absorption spectra according to the highest peak: for phenyl radical, all spectra are shifted by $-437$ cm*-1*, whereas for benzene we introduce different shifts for the on-the-fly ($3010$ cm*-1*), adiabatic harmonic ($3020$ cm*-1*), and vertical harmonic spectra ($3300$ cm*-1*).

IV Results and discussion

IV.1 Absorption spectrum of the phenyl radical

The absorption spectrum of the $\tilde{\text{A}}^{2}\text{B}_{1}\leftarrow\tilde{\text{X}}^{2}\text{A}_{1}$ electronic transition of the phenyl radical has been a subject of theoretical investigation Kim, Mebel, and Lin (2002); Biczysko, Bloino, and Barone (2009); Baiardi, Bloino, and Barone (2013); Patoz, Begušić, and Vaníček (2018) due to a rich vibronic structure originating from the differences between the ground- and excited-state potential energy surfaces. Whereas the adiabatic harmonic approach, in which the potential is expanded about the excited-state minimum, describes the main features of the absorption spectrum, the vertical harmonic model, which expands the potential about the ground-state minimum, fails completely, as shown in our previous work Patoz, Begušić, and Vaníček (2018). Here we compute the phenyl radical absorption spectrum with the on-the-fly ab initio 3TGA with different values of the displacement parameter $d$ and compare the results to the extended TGA spectrum. Since the anharmonicity of the excited-state potential of the phenyl radical is shown to influence the spectrum more than the Herzberg–Teller contribution Patoz, Begušić, and Vaníček (2018), this system is suited for further investigation using the 3TGA, which is specifically intended for treating the anharmonicity of the Herzberg–Teller active modes. The Herzberg–Teller component of the wavepacket has a larger spread in position than does the Condon component, and therefore is more likely to be influenced strongly by the anharmonicity of the potential. Nevertheless, the calculated 3TGA spectra (see Fig. 2) based on $7$ independent trajectories overlap almost perfectly with the ETGA result based on a single trajectory.

This is somewhat surprising considering that several Herzberg–Teller modes are also displaced and, therefore, experience significant anharmonicity of the potential. The results imply that the local harmonic approximation around the true center of the wavepacket is valid despite the additional spread and the change in the shape of the wavepacket. Interestingly, the spectra evaluated with the 3TGA show almost no dependence on the choice of the initial displacement parameter $d$. This is encouraging since a strong dependence would render the 3TGA method impractical; the fact that the dependence is not only weak but almost nonexistent confirms the results of the single-trajectory ETGA in the phenyl radical.

IV.2 Absorption spectrum of benzene

The symmetry-forbidden $\text{B}_{2\text{u}}\leftarrow\text{A}_{1\text{g}}$ transition of benzene provides a beautiful example of a Herzberg–Teller spectrum that is mentioned in many textbooks due to its simple qualitative interpretation Herzberg (1966); Hollas (2004); Quack and Merkt (2011); Bernath (2005). The experimental and theoretical work on the first excited electronic state of benzene Sobolewski and Domcke (1991); Sobolewski, Woywod, and Domcke (1993) and on the vibronic structure of the corresponding electronic spectrum is so extensive that we can only refer to a small fraction here Sponer et al. (1939); Atkinson and Parmenter (1978); Fischer and Jakobson (1979); Faulkner and Richardson (1979); Fischer and Knight (1992); Trost, Stutz, and Platt (1997); Etzkorn et al. (1999); Bernhardsson et al. (2000); He and Pollak (2001); Borges et al. (2003); Schmied et al. (2004); Worth (2007); Loginov, Braun, and Drabbels (2008); Fally, Carleer, and Vandaele (2009); Penfold and Worth (2009); Li et al. (2010); Crespo-Otero and Barbatti (2012). The main progression corresponds to a totally symmetric ring-breathing vibration Herzberg (1966). Although this transition is forbidden by symmetry, its observation in the absorption spectrum is attributed to the non-totally symmetric vibrations which transform as the $e_{2g}$ irreducible representation. Li et al. Li et al. (2010) include the undisplaced distorted modes in the calculation of the absorption spectrum using a global harmonic model without Duschinsky rotation, i.e., by assuming that the ground- and excited-state normal modes are the same. To go beyond the global harmonic model, Penfold and Worth Worth (2007); Penfold and Worth (2009) combine the construction of an anharmonic potential with the multiconfigurational time-dependent Hartree wavepacket dynamics, producing excellent results, but at the cost of a rather detailed inspection of the system required for reducing the computational cost.

Recently, the on-the-fly ab initio ETGA method has been applied to evaluate the absorption spectrum of benzene with rather high accuracy Patoz, Begušić, and Vaníček (2018). As explained above, this method is an automated, simple to use single-trajectory method, and the fact that the ETGA was much more accurate than both vertical and adiabatic global harmonic models is encouraging. Analogous results for the 3TGA are, therefore, presented in Fig. 3, which compares the on-the-fly ab initio result with the global harmonic spectra (top panel). Note that the failure of the vertical harmonic model in the phenyl radical and benzene is not a rule—there are many cases, such as the absorption and photoelectron spectra of ammonia Domcke et al. (1977); Wehrle, Šulc, and Vaníček (2014); Wehrle, Oberli, and Vaníček (2015), in which the vertical is more accurate than the adiabatic harmonic model. Indeed, the vertical harmonic model describes the Franck-Condon region better and so might be expected to be a good model for vertical transitions. The on-the-fly ab initio approach, however, overcomes the issue of choosing the geometry for expanding the potential and is expected to be always at least as good as the better of the two global harmonic models. Moreover, as Fig. 3 shows, the inclusion of anharmonicity in benzene is essential for obtaining a quantitative agreement—compare the 3TGA with the adiabatic harmonic model, which is at least qualitatively correct here.

Finally, the bottom panel of Fig. 3 demonstrates that in benzene the Herzberg–Teller term, captured with the 3TGA, is responsible for the observation of the spectrum because the Condon approximation yields a zero spectrum. This is, indeed, one of the main reasons why we have implemented the 3TGA; the original TGA is often sufficient for Condon spectra.

To investigate the importance of anharmonicity in more detail, the spectra calculated with the 3TGA and ETGA are compared in Fig. 4. Whereas the 3TGA method with the smaller displacement parameter $d$ gives almost the same result as the ETGA approach, the 3TGA spectrum calculated with $d=1$ is red-shifted by $\approx 10$ cm*-1* and the intensity of the first stronger peak is slightly improved. Due to symmetry, the Herzberg–Teller modes cannot be displaced from the minimum, so they do not contribute significantly to the shape of the spectrum. Nevertheless, these modes can influence the total energy shift of the spectrum, which is in accord with the observed results.

IV.3 Norm conservation and fidelity

In contrast with the single-trajectory ETGA, the norm of the 3TGA wavepacket is not conserved because the 3TGA wavepacket is a superposition of wavepackets, each of which feels its own time-dependent potential. Indeed, the numerical results for the phenyl radical (see Fig. 5) and benzene (see Fig. 6) confirm this theoretical prediction from Subsection II.6. Interestingly, larger displacement ($d=1$) , which corresponds to placing the additional Gaussians at the extrema of the initial wavepacket, gives smaller deviations of the norm from unity. In Eq. (42) the time-dependent terms are multiplied by the scaling factor $f_{d}$, indicating that the norm should vary more for smaller displacements, since the factor $f_{d}$ decreases with the displacement. This is in contrast with the expectation that the wavepacket propagation, including the conservation of the norm, should converge to the ETGA with smaller displacements. The fidelity between the 3TGA and ETGA wavepackets,

[TABLE]

can be used as a measure for comparing the quantum dynamics obtained with the two semiclassical approximations; it shows deviations of the three thawed Gaussians from the extended thawed Gaussian wavepacket propagated on the excited-state potential energy surface. As can be seen in Figs. 5 and 6, the fidelity decays similarly for all values of the displacement parameter, despite the differences in the initial overlaps. Although the fidelity decays quickly over time, the final spectra evaluated with the two approximations are nearly the same.

V Conclusion

We have presented the 3TGA, constructed by replacing the Herzberg–Teller part of the initial wavepacket with two displaced Gaussians, in order to describe electronic spectra beyond the Condon approximation and anharmonicity effects beyond the single-trajectory ETGA.

The 3TGA presented here is still a very rough approximation to the exact nuclear wavepacket dynamics. Nevertheless, compared to the global harmonic models, which are frequently employed in computational chemistry, the on-the-fly ab initio methods based on the TGA offer a rather computationally cheap way to partially include anharmonicity. Moreover, the 3TGA allows us to analyze the validity of the single-trajectory ETGA; our results confirm that the additional spread of the wavepacket, induced by the Herzberg–Teller contribution, does not lead to a significant wavepacket splitting during the dynamics induced by electronic absorption in the phenyl radical and benzene.

To conclude, the on-the-fly ab initio implementation of the 3TGA is capable of describing anharmonicity effects and Herzberg–Teller contribution to the spectra, as well as of evaluating both symmetry-forbidden and weakly allowed spectra. The extension beyond the single-trajectory ETGA approach suggests the viability of simple semiclassical methods that can include wavepacket splitting and associated interference, and which would readily outperform global harmonic and single-trajectory methods, while remaining computationally feasible compared to the more advanced semiclassical and quantum methods. Indeed, a related pragmatic approach of using a small number of well-defined, rather than sampled initial conditions, was employed in the multiple coherent states time-averaged semiclassical initial value representation (MC-TA-SC-IVR) Gabas, Conte, and Ceotto (2017) and its “divide-and-conquer” extension Ceotto, Di Liberto, and Conte (2017), which were used to evaluate positions of peaks in vibrational spectra with rather high accuracy. Another appealing feature of the 3TGA is that, despite its simplicity, the resulting spectrum contains the information not only about the positions but also about the intensities of the peaks; this is because the initial conditions are determined by the shape of the initial wavepacket. Last but not least, like MC-TA-SC-IVR, but in contrast to some semiclassical methods that require tens of thousands of trajectories for convergence, the 3TGA, by using only three uniquely defined trajectories (we recommend using $d=1$ always) does not destroy the appealing intuitive picture provided by the original or extended TGA based on a single trajectory. As a result, we anticipate further development of efficient methods using only a few trajectories for semiclassical propagation.

Acknowledgments

The authors acknowledge the financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 683069 – MOLEQULE), from the Swiss National Science Foundation through the NCCR MUST (Molecular Ultrafast Science and Technology) Network, and from the COST action MOLIM (Molecules in Motion). The authors are grateful to Frank Grossmann, David Tannor, and Marius Wehrle for discussions.

Appendix A Orientational average

Let $\overleftrightarrow{T}$ be a second-order tensor in $3$ dimensions (i.e., a $3\times 3$ matrix) and let us evaluate the average of the scalar

[TABLE]

over all unit vectors $\vec{\epsilon}$, i.e.,

[TABLE]

($\int_{S^{2}}\cdots d\Omega$ denotes the integration over the two-dimensional unit sphere $S^{2}$.) We shall prove that the result of this average is simply one third of the trace of the tensor:

[TABLE]

In particular, in Cartesian coordinates,

[TABLE]

Proof: In polar coordinates $(\theta,\varphi)$, the unit vector $\vec{\epsilon}$ is expressed as

[TABLE]

Using Eq. (48), scalar (44) becomes

[TABLE]

where $\cdots$ denote additional, analogous cross terms for $T_{xz}+T_{zx}$ and $T_{yz}+T_{zy}$. The integration over the unit sphere,

[TABLE]

suppresses all the cross terms; e.g., $\sin\varphi\cos\varphi=(1/2)\sin\left(2\varphi\right)$ and the integration over $\varphi$ of the cross term for $T_{xy}+T_{yx}$ is, therefore, zero. Integration of the diagonal terms over $\varphi$ and making a substitution $u:=\cos\theta$ gives

[TABLE]

where in the last step we used the invariance of the trace of a tensor under coordinate transformations, in order to obtain the coordinate-independent expression (46).

Appendix B Proof of the expression (11)

for the gradient of the transition dipole moment

Before it is expressed in the basis of electronic states, the molecular electric dipole operator is simply a sum

[TABLE]

of its nuclear component $\vec{\mu}_{\text{nu}}(q)$ given in Eq. (14) and electronic component

[TABLE]

where $\vec{r}_{n}$ are the Cartesian coordinates of the $n$th electron and $\{\vec{r}_{n}\}:=(\vec{r}_{1},\ldots,\vec{r}_{N_{\text{el}}})$. Relation (11) is proven directly by taking the gradient of Eq. (12):

[TABLE]

The first equality follows from the Leibniz law; the second equality employs a resolution of identity over electronic states and takes into account that the electronic dipole moment operator (49) is independent of nuclear coordinates; the third equality follows from the definition (13) of nonadiabatic coupling vectors, the fact that $\partial_{j}\vec{\mu}_{\text{nu}}(q)$ is a purely nuclear operator, and orthogonality of the electronic states; the fourth equality uses the antihermitian property of the matrix of nonadiabatic couplings,

[TABLE]

which follows from the orthogonality of the electronic states,

[TABLE]

the fifth step completes the proof.

Appendix C Propagation of the width and phase of the

thawed Gaussian wavepacket

The approach by Lee and Heller Lee and Heller (1982) suggests splitting the complex symmetric width matrix $A$ into a product of two matrices $P$ and $Z$:

[TABLE]

By imposing $\dot{Z}_{t}=m^{-1}\cdot P_{t}$, with the initial conditions $Z_{0}=I$, where $I$ is the identity matrix, and $P_{0}=2i\hbar A_{0}$, the expression for the propagation of $Z$ and $P$ parameters is obtained:

[TABLE]

Regarding $\gamma_{t}$, which is a generalization of classical action and represents both the dynamical phase and normalization, it is evaluated as

[TABLE]

where the conditions imposed on $Z$ and $P$ are used in order to obtain the final expression. Note that, since the determinant in the final expression is complex, a proper branch of the logarithm has to be taken in order to make $\gamma_{t}$ continuous in time. If continuity were not imposed on $\gamma_{t}$, the wavepacket would show sudden jumps by $\pi$ in the overall phase. Phase continuity is also important in the evaluation of the correlation function, which comprises a square root of a complex determinant $\det(A_{0}+A_{t}^{\ast})$. The continuity of the correlation function is enforced by taking the appropriate branch of the square root.

Appendix D Derivation of the extended thawed Gaussian

approximation

**Proof of Eq. (29): **The effective potential (21) depends on $t$ only implicitly, via $q_{t}$. Considering $V_{\text{eff}}$ as a function of $q$ and $q_{t}$, its dependence on $q_{t}$ is

[TABLE]

because in the derivation of Eq. (60) the sum $\partial_{q}V|_{q_{t}}+\operatorname{Hess}_{q}V|_{q_{t}}\cdot(q-q_{t})$ appears twice, with opposite signs, and therefore cancels. [Here $V^{\prime\prime\prime}$ is a rank-$3$ tensor of the third derivatives of $V$ and the dot $\cdot$ in the third argument of $V^{\prime\prime\prime}$ indicates that the tensor has been only partially contracted; the right-hand side of Eq. (60) is still a vector.] Within local harmonic approximation (LHA), where all derivatives of $V$ at $q_{t}$ beyond the second are neglected,

[TABLE]

Equation (61) implies that the effective Hamiltonian operator (20), considered as a function of initial conditions $q_{0}$ and $p_{0}$, $\hat{H}_{\text{eff}}(q_{0},p_{0},t)\equiv H_{\text{eff}}(\hat{q},\hat{p},q_{0},p_{0},t)$, is independent of $p_{0}$:

[TABLE]

The effective time evolution operator

[TABLE]

induced by $\hat{H}_{\text{eff}}$ is, in general, a function of $q_{0}$ and $p_{0}$, but factorizing the time-ordered product into a product of exponentials for infinitesimal time steps, expanding the exponential for each time step into a Taylor series, and using Eq. (62) shows that, within the LHA,

[TABLE]

and, indeed, any polynomial in derivatives with respect to $p_{0}$ acting on $\hat{U}$ vanishes:

[TABLE]

The initial Herzberg-Teller state $|\phi_{0}\rangle$ propagated to time $t$ is

[TABLE]

where the Dirac kets $|\psi_{0}\rangle$ and $|\psi_{t}\rangle$ are considered as functions of $q_{0}$ and $p_{0},$ and Eq. (64) was used to switch the order of $P$ and $\hat{U}$. This completes the proof of Eq. (29).

**Proof of Eqs. (30) and (31): **Let us evaluate the derivative $\partial|\psi_{t}\rangle/\partial p_{0}$, needed in Eq. (65), analytically in position representation, in which $|\psi_{t}\rangle$ is the TGA wavefunction (19):

[TABLE]

In the first step of the derivation, we neglected the dependence of the width matrix $A_{t}$ on $p_{0}$, which follows from Eq. (24) within the LHA. In the second step, we used the derivative

[TABLE]

where $S_{t}=\int_{0}^{t}L_{\tau}d\tau$ is the classical action component of $\gamma_{t}$ expressed either as a function of $q_{0}$ and $p_{0}$ or of $q_{0}$ and $q_{t}$. The derivation of the central ETGA equations (30) and (31) is completed by substituting the Herzberg-Teller form [Eq. (10)] of the polynomial $P$,

[TABLE]

into Eq. (29), giving

[TABLE]

and using expression (66) for $\partial\psi_{t}(q)/\partial p_{0}$.

Appendix E Extrema of the Herzberg–Teller part of

the wavepacket

The extrema of the wavepacket (32) are found by setting its gradient to zero. The gradient, given by

[TABLE]

where $\Delta q=q-q_{0}$, will be zero if and only if

[TABLE]

because the Gaussian function $g_{q_{0}}(q)$ is strictly positive. Since the initial width matrix $A_{0}$ is positive definite, it has an inverse and we can multiply Eq. (68) on the left with $A_{0}^{-1}$, which gives

[TABLE]

Scalar product of Eq. (69) with the vector $\partial_{q}\mu$ gives

[TABLE]

with two solutions

[TABLE]

Finally, substitution of $\partial_{q}\mu^{T}\cdot\Delta q$ from Eq. (71) into Eq. (69) yields

[TABLE]

where it is easy to see that the positive sign corresponds to the local maximum and negative sign to the local minimum, completing the proof of Eqs. (35) and (36).