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

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.


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 [1,2], 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 [3,4], 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 [5,6].
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 [7][8][9][10][11][12][13][14][15], one of the simplest, yet very intuitive and often surprisingly accurate, is the thawed Gaussian approximation (TGA) of Heller [16]. Although his work on Gaussian wavepackets [17] inspired a number of "direct", i.e., on-the-fly ab initio quantum [18][19][20][21][22] and semiclassical [23][24][25][26][27][28] methods, the original thawed Gaussian approximation was only revisited recently, when it was implemented within the onthe-fly ab initio framework for evaluating molecular absorption, emission, and photoelectron spectra [29,30]. Due to its computational efficiency in comparison with exact quantum dynamics methods, the TGA can account for [29,30] the full dimensionality of the system and mode-mixing (Duschinsky effect), both of which can significantly influence the wavepacket dynamics [31]. Recently, we have implemented [32] 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) [32,33], propagates a Gaussian wavepacket multiplied by a general polynomial in nuclear coordinates using a local harmonic approximation, i.e., the secondorder 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 R D , electronic S-dimensional complex Hilbert space C S , and the ambient 3-dimensional space R 3 .
To distinguish these spaces, 3-dimensional vectors will be denoted with an arrow (e.g., ), whereas D-dimensional nuclear vectors or D × 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 µ 21 · or p T · m −1 · p).
We shall use the bold font for S × S matrices (such as µ) 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 AB).

II. THEORY
A. 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 ω can be expressed as the Fourier transform of the dipole time autocorrelation function C µµ ( , t). Assuming the zero temperature approximation, the initial state is |1, g , 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 [3] C µµ ( , t) = 1, g|e iĤ 1 t/ μ 12 e −iĤ 2 t/ μ 21 |1, g , whereĤ 1 andĤ 2 are nuclear Hamiltonian operators in the ground and excited electronic states, andμ is the projection of the matrix elementˆ µ 21 of the molecular transition dipole moment matrix µ on a three-dimensional polarization unit vector of the electric field.

B. 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 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 [34,35].
It is useful to define the spectrum tensor ← → σ (ω), from which the spectrum (1) for a specific polarization , is obtained by "evaluation": The orientational averaging of the spectrum corresponds to the averaging of σ( , ω) over all unit vectors . 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: Appendix A contains an explicit proof of these equalities, of which the first is coordinateindependent and the other two are explicit in Cartesian coordinates ( 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: Now let us go back to our two-state system, where only one element of the electric transition dipole moment, µ 21 (q) plays a role and C µµ ( , t) is given by Eq. (2). Then it is convenient to suppress the subscripts 21 and denote the 21 element simply as µ(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 µµ ( , t) either for three orthogonal polarizations and a fixed molecular orientation [36], or for three orthogonal orientations of the molecule and a fixed polarization . Within the Condon approximation, in which µ is coordinate-independent, the orientational average (7) simplifies further into the standard textbook recipe where | µ| is the magnitude of the transition dipole moment, so only a single calculation is required-for the transition dipole moment aligned with the field.

C. 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 [37], this moment is assumed to be independent of the molecular geometry and is approximated to the zeroth order around the initial geometry: 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 α ← β with µ αβ (q 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 [38] that takes into account also the gradient of the transition dipole moment with respect to nuclear degrees of freedom: where ∂ q is an abbreviation for ∂/∂q. Although the Herzberg-Teller approximation is often discussed in terms of vibronic couplings between different electronic states [39], this relation is not obvious from Eq. (10).
Here we demonstrate this relationship explicitly. In particular, we shall prove a remarkable equality ∂ j µ(q) = [ µ(q), F j (q)] + ∂ j µ nu (q)1 (11) satisfied by the gradient of the matrix representation µ(q) of the molecular dipole operator µ mol at the nuclear configuration q. In this equation, matrix elements of µ(q) are defined as partial, electronic expectation values elements of the matrix F j (q) of nonadiabatic vector couplings are obtained as and µ nu is the nuclear component ofˆ µ mol : Here Z N is the atomic number and 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 con- that explains the meaning of intensity borrowing, in which the symmetry forbidden transition to the formally dark state α occurs due to a borrowing of the transition intensity from the neighboring bright electronic states γ that are vibronically coupled to the state α.
Note that, despite introducing nonadiabatic couplings between excited electronic states, the original Born-Oppenheimer picture adopted in Subsection II A 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 fieldfree 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.

E. 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 in terms of the vibrational zero-point energy E 1,g of the ground electronic state and the wavepacket autocorrelation function for the initial wavepacket |φ(0) =μ 21 |1, g propagated on the excited-state potential energy surface with the HamiltonianĤ 2 : 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 [16], 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 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 γ t a complex number whose real part is a dynamical phase and imaginary part guarantees the normalization for times where m is the diagonal mass matrix and V eff an effective time-dependent potential given by the local harmonic approximation of the true potential V at the center of the wavepacket: To avoid the confusion between the Laplacian and Hessian operators, both of which are often denoted by ∂ 2 /∂q 2 , we use the notation Hess q V for Hessian matrix with respect to nuclear coordinates: Insertion of the wavepacket ansatz (19) and the effective potential (21) into the timedependent Schrödinger equation yields a system of ordinary differential equations [16] 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 γ t involves propagating the 2D × 2D stability matrix 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 [33] proposed the extended TGA (ETGA) that considers a more general form of the initial wavepacket, namely a Gaussian multiplied by a polynomial [32,33]: 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 ψ 0 (q) on p 0 comes from the exponent p T 0 · (q − q 0 ) [see Eq. (19)], the polynomial prefactor in Eq. (27) can be replaced by the same polynomial in the derivatives with respect to p 0 : In Appendix D, we prove that the local harmonic approximation implies that the ETGA wavepacket at time t has the simple form [33] φ t (q) = P i where the dependence of ψ 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 γ 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 where the linear parameter of the polynomial at time t is

F. 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 µ 21 (q 0 ) = 0], where the probability density at the center is zero. Although the ETGA is still an exact solution of the timedependent Schrödinger equation in a global harmonic potential, its performance in a general potential can, therefore, be questioned [33]. To further investigate the errors introduced by the local harmonic approximation, we propose to approximate the Herzberg-Teller part of the initial state used in the ETGA with an antisymmetric linear combination of two displaced Gaussians (see Fig. 1): where g qc (q) is a normalized Gaussian ψ 0 (q) centered at q c instead of q 0 , with zero initial momentum and phase (p 0 = γ 0 = 0): ∆ d is a displacement vector and f d a scaling factor. Note also that the same width matrix A 0 is used for g qc as for ψ 0 (q).
We will now show that the approximative wavepacket φ 3TGA-HT 0 converges to the exact Herzberg-Teller wavepacket φ ETGA-HT 0 when the displacement ∆ 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 µ instead of First, let us discuss the shape of the exact wavepacket φ ETGA-HT 0 . Obviously, in one nuclear dimension, φ ETGA-HT 0 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 with the displacement vector It is rather obvious that in order for φ 3TGA-HT 0 to be a good fit of φ ETGA-HT 0 , the two Gaussians in Eq. (33) must be displaced from q 0 along the direction of the extrema of φ ETGA-HT 0 . We therefore choose where the dimensionless parameter d controls the magnitude of the displacement. The , which gives A quick inspection reveals φ ETGA-HT 0 to be a directional derivative of g q 0 (q) and φ 3TGA-HT 0 a finite-difference approximation of this derivative, exact in the limit d → 0. One can see this explicitly by noting that the normalized initial overlap of the two wavepackets, 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 To see this, consider several wavepackets |ψ i , each propagated with its own evolution op-eratorÛ i . IfÛ i =Û j , then the overlap of two such wavepackets is time-dependent: The time dependence of the norm of the 3TGA wavepacket arises from these time-dependent overlaps since Due to the time dependence of the norm, the calculated dipole time autocorrelation function has to be renormalized at each time step.

G. Ab initio implementation
The ab initio implementation of the thawed Gaussian approximation has been discussed in Refs. 29 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 (π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
A. Absorption spectrum of the phenyl radical The absorption spectrum of theÃ 2 B 1 ←X 2 A 1 electronic transition of the phenyl radical has been a subject of theoretical investigation [32,[42][43][44] 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 [32]. 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 [32], 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.
The main progression corresponds to a totally symmetric ring-breathing vibration [2]. 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. [40] 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 [61,64] 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 [32]. 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 [29,30,66], 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 ≈ 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.

C. 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 F. ing 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, 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) [28] and its "divide-and-conquer" extension [14], 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.
( S 2 · · · dΩ 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: In particular, in Cartesian coordinates, Proof : In polar coordinates (θ, ϕ), the unit vector is expressed as = (sin θ cos ϕ, sin θ sin ϕ, cos θ) .
Using Eq. (A5), scalar (A1) becomes T ( ) = T xx cos 2 ϕ + T yy sin 2 ϕ sin 2 θ + T zz cos 2 θ where · · · denote additional, analogous cross terms for T xz + T zx and T yz + T zy . The integration over the unit sphere, suppresses all the cross terms; e.g., sin ϕ cos ϕ = (1/2) sin (2ϕ) and the integration over ϕ of the cross term for T xy + T yx is, therefore, zero. Integration of the diagonal terms over ϕ and making a substitution u := cos θ gives 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 (A3).
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ˆ of its nuclear component µ nu (q) given in Eq. (14) and electronic component where r n are the Cartesian coordinates of the nth electron and { r n } := ( r 1 , . . . , r N el ). Relation (11) is proven directly by taking the gradient of Eq. (12): 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 (B1) is independent of nuclear coordinates; the third equality follows from the definition (13) of nonadiabatic coupling vectors, the fact that ∂ j µ 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, which follows from the orthogonality of the electronic states, 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 [33] suggests splitting the complex symmetric width matrix A into a product of two matrices P and Z: By imposingŻ t = m −1 ·P t , with the initial conditions Z 0 = I, where I is the identity matrix, and P 0 = 2i A 0 , the expression for the propagation of Z and P parameters is obtained: Regarding γ t , which is a generalization of classical action and represents both the dynamical phase and normalization, it is evaluated as 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 γ t continuous in time. If continuity were not imposed on γ t , the wavepacket would show sudden jumps by π 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 ). The continuity of the correlation function is enforced by taking the appropriate branch of the square root. Considering V eff as a function of q and q t , its dependence on q t is because in the derivation of Eq. (D1) the sum ∂ q V | qt +Hess q V | qt ·(q −q t ) appears twice, with opposite signs, and therefore cancels. The initial Herzberg-Teller state |φ 0 propagated to time t is |φ t =Û (q 0 , p 0 , t)|φ 0 =Û (q 0 , p 0 , t)P i ∂ ∂p 0 |ψ 0 = P i ∂ ∂p 0 Û (q 0 , p 0 , t)|ψ 0 = P i where the Dirac kets |ψ 0 and |ψ t are considered as functions of q 0 and p 0 , and Eq. (D5) was used to switch the order of P andÛ . This completes the proof of Eq. (29).
Proof of Eqs. (30) and (31): Let us evaluate the derivative ∂|ψ t /∂p 0 , needed in Eq. (D6), analytically in position representation, in which |ψ t is the TGA wavefunction (19): 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 where S t = t 0 L τ dτ is the classical action component of γ 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 , into Eq. (29), giving φ t (q) = µ 21 (q 0 ) + ∂ q µ 21 | T q 0 · i ∂ ∂p 0 ψ t (q), and using expression (D7) for ∂ψ t (q)/∂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 where ∆q = q − q 0 , will be zero if and only if ∂ q µ = 2(∂ q µ T · ∆q)A 0 · ∆q (E2) 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. (E2) on the left with A −1 0 , which gives A −1 0 · ∂ q µ = 2(∂ q µ T · ∆q)∆q (E3) Scalar product of Eq. (E3) with the vector ∂ q µ gives with two solutions Finally, substitution of ∂ q µ T · ∆q from Eq. (E5) into Eq. (E3) yields ∆q = ± √ 2 2