Single-Hessian thawed Gaussian approximation: The missing rung on the ladder

TL;DR

This paper introduces a computationally efficient single-Hessian thawed Gaussian approximation for molecular spectra that maintains energy conservation and performs comparably to more expensive methods, improving semiclassical simulations.

Contribution

The paper presents a novel single-Hessian thawed Gaussian method that conserves energy and approximates potential Hessians with a constant matrix, enhancing efficiency in semiclassical calculations.

Findings

Performs nearly as well as on-the-fly ab initio methods

Significantly better than global harmonic schemes

Conserves energy exactly despite time-dependent Hamiltonian

Abstract

To alleviate the computational cost associated with on-the-fly ab initio semiclassical calculations of molecular spectra, we propose the single-Hessian thawed Gaussian approximation, in which the Hessian of the potential energy at all points along an anharmonic classical trajectory is approximated by a constant matrix. The spectra obtained with this approximation are compared with the exact quantum spectra of a one-dimensional Morse potential and with the experimental spectra of ammonia and quinquethiophene. In all cases, the single-Hessian version performs almost as well as the much more expensive on-the-fly ab initio thawed Gaussian approximation and significantly better than the global harmonic schemes. Remarkably, unlike the thawed Gaussian approximation, the proposed method conserves energy exactly, despite the time dependence of the corresponding effective Hamiltonian, and, in addition, can be mapped to a higher-dimensional time-independent classical Hamiltonian system. We also provide a detailed comparison with several related approximations used for accelerating prefactor calculations in semiclassical simulations.

Tables

Table 1. Table 1: Equations of motion for parameters A t subscript 𝐴 𝑡 A_{t} and γ t subscript 𝛾 𝑡 \gamma_{t} of the Gaussian wavepacket ( 7 ) propagated with the single-Hessian thawed Gaussian approximation, exact for the approximate Hamiltonian ( 23 ), and with the single-Hessian and single-trajectory versions of the approximations discussed in Appendix  A . α ~ 0 subscript ~ 𝛼 0 \tilde{\alpha}_{0} is defined in Eq. ( 48 ), A ref subscript 𝐴 ref A_{\text{ref}} in Eq. ( 49 ), and we use − ( i ​ ℏ / 2 ) ​ Tr ​ ( m − 1 ⋅ A ref ) = ∑ j = 1 D 1 2 ​ ℏ ​ ω ref , j 𝑖 Planck-constant-over-2-pi 2 Tr ⋅ superscript 𝑚 1 subscript 𝐴 ref superscript subscript 𝑗 1 𝐷 1 2 Planck-constant-over-2-pi subscript 𝜔 ref 𝑗 -(i\hbar/2)\text{Tr}(m^{-1}\cdot A_{\text{ref}})=\sum_{j=1}^{D}\frac{1}{2}\hbar\omega_{\text{ref},j} . For the single-Hessian Herman–Kluk approximation, we assume that matrices Γ Γ \Gamma , m 𝑚 m , and V ref ′′ ​ ( q ref ) subscript superscript 𝑉 ′′ ref subscript 𝑞 ref V^{\prime\prime}_{\text{ref}}(q_{\text{ref}}) commute. For the frozen Gaussian approximation, Heller ( 1981b ); Davis and Heller ( 1981 ) the general expression γ ˙ t = p t T ⋅ m − 1 ⋅ p t − ⟨ H ^ ⟩ subscript ˙ 𝛾 𝑡 ⋅ superscript subscript 𝑝 𝑡 𝑇 superscript 𝑚 1 subscript 𝑝 𝑡 delimited-⟨⟩ ^ 𝐻 \dot{\gamma}_{t}=p_{t}^{T}\cdot m^{-1}\cdot p_{t}-\langle\hat{H}\rangle is expanded using the single-Hessian potential ( 23 ) and the total energy of a Gaussian wavepacket [Eqs. ( 59 )–( 61 )] applied to a coherent state ( A t = A 0 subscript 𝐴 𝑡 subscript 𝐴 0 A_{t}=A_{0} ).

Thawed Gaussian	Frozen Gaussian
Herman–Kluk	Adiabatic Herman–Kluk
Johnson	Prefactor-free

Table 2. Table 2: Cosines of the spectral contrast angles [Eq. ( 39 )] comparing the exact spectrum of a Morse potential ( χ = 0.02 𝜒 0.02 \chi=0.02 , see Sec.  III.1 ) with the spectra evaluated using the single-Hessian thawed Gaussian approximation (SH TGA) and single-Hessian single-trajectory approximations discussed in Appendix  A . In a one-dimensional system, the adiabatic Herman–Kluk approximation is equivalent to the Herman–Kluk method. Results for adiabatic, vertical, and initial reference Hessians are shown. The top three rows contain contrast angles of the spectra shifted so that their overlaps with the exact (reference) spectrum are maximal; the rows below refer to the unshifted spectra, where the errors due to constant horizontal shifts of the spectra are accounted for.

Reference	SH	Herman	Johnson	Frozen	Prefactor
Hessian	TGA	Kluk		Gaussian	free
			Shifted
Adiabatic	0.975	0.973	0.973	0.973	0.973
Vertical	0.964	0.973	0.973	0.973	0.973
Initial	0.973	0.973	0.973	0.973	0.973
			Not shifted
Adiabatic	0.975	0.973	0.973	0.973	0.006
Vertical	0.242	0.243	0.243	0.172	0.006
Initial	0.899	0.899	0.899	0.899	0.006

Equations

$$|\psi (t)⟩=e^{-i\hat{H}t/\hslash }|{\psi }_i⟩,$$

$$\hat{H}=H(\hat{q},\hat{p})=\frac{1}{2}{\hat{p}}^T\cdot m^{-1}\cdot \hat{p}+V(\hat{q}),$$

$$C(t)=⟨{\psi }_i|\psi (t)⟩.$$

$${\sigma }_{\text{abs}}(\omega )=\frac{4\pi \omega }{3\hslash c}|{\vec{\mu }}_{21}{|}^2\text{Re}{\int }_0^{\infty }C(t)e^{i(\omega +E_{1,i}/\hslash )t}dt,$$

$${\sigma }_{\text{em}}(\omega )=\frac{4{\omega }^3}{3\pi \hslash c^3}|{\vec{\mu }}_{21}{|}^2\text{Re}{\int }_0^{\infty }C(t{)}^{\ast }{e}^{i(\omega -E_{2,i}/\hslash )t}dt,$$

$${\sigma }_{\text{abs}}(\omega )=\frac{4{\pi }^2\omega }{3\hslash c}|{\vec{\mu }}_{21}{|}^2\sum _n|⟨n|{\psi }_i⟩{|}^2\delta (\omega -\frac{E_{2,n}-E_{1,i}}{\hslash }),$$

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

$$V_{\text{LHA}}(q,t)=V(q_t)+V^{\prime }(q_t{)}^T\cdot (q-q_t)+\frac{1}{2}(q-q_t{)}^T\cdot V^{\prime \prime }(q_t)\cdot (q-q_t),$$

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

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

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

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

$$\frac{dE}{dt}$$

$$=⟨\psi (t)|\frac{d}{dt}{\hat{H}}_{\text{eff}}(t)|\psi (t)⟩$$

$$=⟨\psi (t)|\frac{d}{dt}{\hat{V}}_{\text{LHA}}(t)|\psi (t)⟩$$

$$=\frac{1}{2}⟨\psi (t)|(\hat{q}-q_t{)}^T\cdot B_t\cdot (\hat{q}-q_t)|\psi (t)⟩$$

$$=\frac{1}{2}\mathrm{Tr}(B_t\cdot {\Sigma }_t^2),$$

$${\Sigma }_t^2$$

$$=\int dq|\psi (q,t){|}^2(q-q_t)\otimes (q-q_t{)}^T$$

$$\displaystyle=\bigg{(}\frac{2}{\hbar}\operatorname{Im}A_{t}\bigg{)}^{-1}$$

$$\frac{d}{dt}=p_t^T\cdot m^{-1}\cdot \frac{d}{d{q}_t}$$

$$V_{\text{HA}}(q)=V_{\text{eq}}+\frac{1}{2}(q-q_{\text{eq}}{)}^T\cdot k\cdot (q-q_{\text{eq}}).$$

$$V_{\text{SH}}(q,t)=V(q_t)+V^{\prime }(q_t{)}^T\cdot (q-q_t)+\frac{1}{2}(q-q_t{)}^T\cdot V_{\text{ref}}^{\prime \prime }(q_{\text{ref}})\cdot (q-q_t).$$

$$\frac{dE}{dt}$$

$$=⟨\psi (t)|b_t^T\cdot (\hat{q}-q_t)|\psi (t)⟩$$

$$=b_t^T\cdot ⟨\psi (t)|\hat{q}-q_t|\psi (t)⟩$$

$$=0,$$

$$A_t$$

$$P_t$$

$${\dot{Q}}_t$$

$${\dot{P}}_t$$

$${\dot{Q}}_t=\frac{\partial H_{\text{sc}}}{\partial P^{\ast }}\text{ and }{\dot{P}}_t=-\frac{\partial H_{\text{sc}}}{\partial Q^{\ast }}$$

$$H_{\text{sc}}(Q,P;q_t)=\frac{1}{2}\mathrm{Tr}[P^{†}\cdot m^{-1}\cdot P+Q^{†}\cdot V^{\prime \prime }(q_t)\cdot Q],$$

$$H_{\text{map}}(q,p,Q,P)=H(q,p)+H_{\text{sc}}(Q,P).$$

$$V(q)=V_{\text{eq}}+D_e[1-e^{-a(q-q_{\text{eq}})}{]}^2,$$

$${\omega }_{\text{eq}}=\sqrt{V^{\prime \prime }(q_{\text{eq}})}=\sqrt{2D_e{a}^2},$$

$$\chi =\frac{{\omega }_{\text{eq}}}{4D_e}.$$

Full text

I Introduction

Simulation of vibrationally resolved electronic spectra of large polyatomic molecules is a challenge for computational chemistry. The exact calculation is impossible for most but smallest molecular systems due to the exponentially scaling cost of computing the full potential energy surfaces of the electronic states involved in the transition. In the well-known time-independent formalism, the intensities of the individual vibronic transitions are determined by the Franck–Condon factors, i.e., the squares of overlaps between the vibrational eigenstates of the two electronic states, while the frequencies of transitions are given by the differences of the corresponding vibrational eigenvalues. A popular method for computing vibronic spectra constructs global harmonic models of the two potential energy surfaces.Santoro et al. (2007, 2008); Barone et al. (2009) Then, the vibrational functions, as well as their overlaps, are given analytically. Anharmonic corrections can be included perturbativelyBonness et al. (2006); Luis, Bishop, and Kirtman (2004); Yang et al. (2012); Egidi et al. (2017) or variationally.Luis, Kirtman, and Christiansen (2006); Bowman, Carrington, and Meyer (2008); Koziol et al. (2009); Meier and Rauhut (2015) In smaller systems, it is feasible to apply anharmonic corrections to both eigenstates and eigenvalues, which affects both positions and intensities of vibronic transitions.Mok et al. (2000); Luis, Bishop, and Kirtman (2004); Koziol et al. (2009) In larger systems, however, this is computationally challenging and the anharmonic corrections are almost exclusively included only through the frequencies, without affecting the Franck–Condon factors.Egidi et al. (2014, 2017); Biczysko, Krupa, and Wierzejewska (2018)

Time-dependent approaches, based on computing the dipole time correlation function,Heller (1981a); Mukamel (1999); Tannor (2007) have also been developed at different levels of accuracy, ranging from global harmonic modelsBaiardi, Bloino, and Barone (2013) to exact quantum dynamics methodsMeyer, Gatti, and Worth (2009) on anharmonic potential energy surfaces. The time-dependent formalism allows for an on-the-fly implementation, where the potential data are evaluated only when needed, and therefore provides an easier route to including anharmonicity. We focus our attention on the thawed Gaussian approximation,Heller (1975); Grossmann (2006) which, as several other semiclassicalTatchen and Pollak (2009); Ceotto et al. (2009a, b); Wong et al. (2011); Gabas, Conte, and Ceotto (2017); Gabas et al. (2018) and quantumBen-Nun, Quenneville, and Martínez (2000); Saita and Shalashilin (2012); Richings et al. (2015) dynamics methods, has been implemented in an on-the-fly fashion and combined with an ab initio evaluation of the potential.Wehrle, Šulc, and Vaníček (2014); Wehrle, Oberli, and Vaníček (2015) The method assumes validity of the Born–Oppenheimer approximation and propagates a Gaussian wavepacket in a locally harmonic potential constructed about the current center of the wavepacket at each time step. This rather simple propagation scheme, proposed by Heller as the first step beyond the global harmonic approximation in the hierarchy of time-dependent methods, was shown to work well for low or medium resolution electronic spectra, where only short-time propagation of the wavepacket is needed.Wehrle, Šulc, and Vaníček (2014); Wehrle, Oberli, and Vaníček (2015); Patoz, Begušić, and Vaníček (2018) To further reduce the computational cost of on-the-fly ab initio calculations, one can employ a Hessian interpolation scheme, in which the Hessians are evaluated only every several steps and interpolated in between.Wehrle, Šulc, and Vaníček (2014)

Here, we propose a new approach, which still uses a fully anharmonic classical trajectory to guide the Gaussian wavepacket but only a single Hessian to propagate the width. Hence, this “single-Hessian thawed Gaussian approximation” further reduces the cost of spectra calculations to that of a single classical trajectory. The method is validated on a Morse potential as well as on full-dimensional on-the-fly ab initio simulations of the absorption spectrum of ammonia and emission spectrum of quinquethiophene. The single-Hessian method performs better than the global harmonic approaches and in some cases even better than the standard thawed Gaussian approximation. Although the effective Hamiltonian associated to the single-Hessian thawed Gaussian approximation is time-dependent, we demonstrate—both analytically and numerically—that the energy is conserved. Finally, we explore the relation between this single-Hessian approach and similar well-known approximations to the prefactor in the semiclassical Herman–Kluk initial value representation.

II Theory

II.1 Time-dependent approach to vibrationally

resolved electronic spectroscopy

Let $|\psi(t)\rangle$ be a wavepacket

[TABLE]

propagated with a time-independent Hamiltonian

[TABLE]

where $|\psi_{i}\rangle$ represents the initial state. Within the electric dipole approximation, first-order perturbation theory, and assuming the Condon approximation, vibrationally resolved electronic spectra can be computed from the wavepacket autocorrelation function

[TABLE]

The type of spectroscopy determines the choice of $\psi_{i}$ and $H$. If $\psi_{i}$ is a vibrational eigenstate of the ground electronic state $1$ and $H$ the excited-state vibrational Hamiltonian, the rotationally averaged absorption cross-section is evaluated as the Fourier transformHeller (1981a); Lami, Petrongolo, and Santoro (2004); Tannor (2007)

[TABLE]

where $E_{1,i}$ is the energy of state $\psi_{i}$ before photon absorption and $\vec{\mu}_{21}$ the transition dipole moment between the ground and excited electronic states evaluated at the ground-state equilibrium geometry. The emission spectrum, measured as the emission rate per unit frequency, is obtained by taking the $\psi_{i}$ to be the vibrational eigenstate of the excited electronic state $2$ and $H$ the ground-state vibrational Hamiltonian:Lami, Petrongolo, and Santoro (2004); Niu et al. (2010)

[TABLE]

where $E_{2,i}$ is the energy of state $\psi_{i}$ before photon emission. Spectra defined in Eqs. (4) and (5) are positive at all frequencies, which can be shown by inserting a resolution of identity in the expression (3) for the autocorrelation function to derive the time-independent expression; e.g., for the absorption spectrum, one obtainsTannor (2007)

[TABLE]

where $|n\rangle$ are the eigenstates of the excited-state vibrational Hamiltonian with energies $E_{2,n}$. Equations (4) and (6) are equivalent for any time-independent Hamiltonian. However, if the true time-independent Hamiltonian is approximated by an effective time-dependent one, for example, through the local harmonic or cubic approximations, negative spectral features may arise.Wehrle, Oberli, and Vaníček (2015)

II.2 Thawed Gaussian approximation

Evaluation of the autocorrelation function (3) requires propagating the vibrational wavepacket; among many quantum and semiclassical methods, one of the simplest is the thawed Gaussian approximation.Heller (1975) A thawed Gaussian wavepacket is described by its time-dependent position $q_{t}$, momentum $p_{t}$, complex symmetric matrix $A_{t}$, and a complex number $\gamma_{t}$:

[TABLE]

where $N_{0}=[\det(\text{Im}A_{0}/\pi\hbar)]^{1/4}$ is a normalization constant. Classical parameters $q_{t}$ and $p_{t}$ are the expectation values of the position and momentum; the imaginary part of matrix $A_{t}$ controls the width of the wavepacket, while its real part introduces a spatial chirp; the real part of $\gamma_{t}$ is a time-dependent phase factor, while its imaginary part ensures normalization at all times. The Gaussian form (7) is exactly preserved under evolution in a harmonic potential, even a time-dependent one. In the thawed Gaussian approximation, the wavepacket (7) is propagated with an effective Hamiltonian $\hat{H}_{\text{eff}}(t)=\hat{T}+\hat{V}_{\text{LHA}}(t)$ given by the sum of the kinetic energy $T$ and the time-dependent local harmonic approximation $V_{\text{LHA}}$ of the true potential $V$ about $q_{t}$:

[TABLE]

with $V^{\prime}(q_{t})$ representing the gradient and $V^{\prime\prime}(q_{t})$ the $D\times D$ Hessian matrix of the potential evaluated at the center of the wavepacket $q_{t}$. Inserting the Gaussian ansatz (7) and effective potential (8) into the time-dependent Schrödinger equation gives the following equations of motion for the wavepacket parameters:Heller (1975)

[TABLE]

where $m$ is the mass matrix and $L_{t}$ the Lagrangian.

If the wavepacket remains localized, the effective locally harmonic potential is a good approximation and the thawed Gaussian propagation is expected to be rather accurate. The approximation accounts partially for anharmonicity by propagating the wavepacket’s center $(q_{t},p_{t})$ classically with the true, anharmonic potential $V(q)$ [Eqs. (9)–(10) are Hamilton’s equations of motion for $H(q,p)$] and by accounting for the changes in its Hessian, which affect the semiclassical parameters $A_{t}$ and $\gamma_{t}$. Another advantage of the thawed Gaussian approximation is its efficiency: it requires propagating four time-dependent parameters, which depend only on the local potential information.

Yet, there are also several drawbacks: First, the Gaussian ansatz (7) cannot describe wavepacket splitting, tunneling, or nonadiabatic effects. In very anharmonic systems, where the exact wavepacket splits and delocalizes quickly, the thawed Gaussian wavepacket behaves unphysically. Thus, the method is limited to short propagation times and low-resolution electronic spectra. Second, because the effective potential (8) is, in general (i.e., for potentials beyond quadratic), time-dependent, the thawed Gaussian approximation does not conserve energy:Wehrle, Oberli, and Vaníček (2015); Rohrdanz and Cina (2006)

[TABLE]

where $B_{t}:=p_{t}^{T}\cdot m^{-1}\cdot V^{\prime\prime\prime}(q_{t})$, $V^{\prime\prime\prime}(q_{t})$ is a rank-3 tensor of third derivatives of the potential with respect to position, and

[TABLE]

is the position covariance matrix. Equation (14) follows because the thawed Gaussian solves exactly the Schrödinger equation with $H_{\text{eff}}$, while Eq. (15) relies on the time independence of the kinetic energy operator. To derive Eq. (16), we used the chain rule

[TABLE]

for the differentiation of the energy, gradient, and Hessian evaluated at position $q_{t}$. As noted already in Section II.1, the time dependence of the effective Hamiltonian also leads to unphysical negative intensities in the spectra.

II.3 Hessian interpolation

To reduce the cost of ab initio Hessian calculations, the on-the-fly ab initio thawed Gaussian approximation is readily combined with an interpolation scheme, where the Hessians are computed only every few steps and the intermediate Hessians are obtained from a second-order polynomial interpolation. Typically, the Hessians need to be computed only every four to eight time steps.Wehrle, Šulc, and Vaníček (2014); Begušić, Roulet, and Vaníček (2018) Since the Hessians are not needed for the propagation of the classical trajectory, additional speed-up is achieved through parallel computation of the Hessians after the full trajectory is known. Note that other Hessian approximations, such as the Hessian update schemesCeotto, Zhuang, and Hase (2013); Zhuang et al. (2013); Ianconescu, Tatchen, and Pollak (2013) and Gaussian process regressionAlborzpour, Tew, and Habershon (2016); Laude et al. (2018) have been developed in the context of ab initio simulations. The considerable cost of multiple Hessian evaluations has also inspired various semiclassical approximations,Di Liberto and Ceotto (2016) including the prefactor-free,Tatchen and Pollak (2009) adiabatic,Guallar, Batista, and Miller (1999, 2000) harmonic,Di Liberto and Ceotto (2016) and “poor person’s”Tatchen et al. (2011) variations of the Herman–Kluk propagator.

II.4 Global harmonic approximation

In computational chemistry, most calculations of vibrationally resolved electronic spectra employ the global harmonic models, where the true potential energy surface is approximated as

[TABLE]

In Eq. (22), $V_{\text{eq}}$ is the potential energy and $q_{\text{eq}}$ the position of the minimum of the harmonic potential with a force constant matrix $k$. In practice, the global harmonic model is constructed from ab initio data evaluated at a single molecular geometry, which makes such calculations feasible for rather large systems. The thawed Gaussian wavepacket (7) is exact in the harmonic potential (22) and can be propagated analytically. Furthermore, because the potential is time-independent, the energy is conserved exactly and the corresponding spectra do not suffer from unphysical negative intensities. However, the method neglects anharmonicity completely and, therefore, is less accurate than the thawed Gaussian approximation.

II.5 Single-Hessian thawed Gaussian approximation

Let us now consider using a single Hessian in the local harmonic approximation (8), e.g., by choosing a reference Hessian $V_{\text{ref}}^{\prime\prime}(q_{\text{ref}})$ and approximating the potential at each point in time as

[TABLE]

The single-Hessian thawed Gaussian approximation, which propagates the wavepacket (7) in the effective potential (23), is, obviously, even more efficient than the original thawed Gaussian approximation; the single-Hessian analogue requires only one Hessian to be evaluated for the whole propagation, i.e., its cost is almost the same as running a single classical trajectory. Because the effective potential (23) is Hermitian, the single-Hessian method conserves the norm of the wavefunction. As for the accuracy, the approximation (23) of the potential still includes anharmonicity partially through the first two terms and thus is more accurate than the global harmonic approximation, but is clearly worse than the local harmonic approximation (8) (see Fig. 1). Yet, the single-Hessian approach also results in several improvements related to spectra calculations:

First, the propagation of matrix $A_{t}$ is now determined exclusively by the reference Hessian and is decoupled from the classical dynamics of $q_{t}$ and $p_{t}$. Therefore, the wavepacket does not spread or contract unphysically in an attempt to describe wavepacket splitting, but rather stays compact at all times, similarly to a squeezed state in a globally harmonic potential. We show in several numerical examples that this feature is preferred in more anharmonic potentials.

Second, the single-Hessian thawed Gaussian approximation conserves energy exactly:

[TABLE]

where $b_{t}^{T}:=p_{t}^{T}\cdot m^{-1}\cdot(V^{\prime\prime}(q_{t})-V_{\text{ref}}^{\prime\prime}(q_{\text{ref}}))$. Above, we used the time independence of the kinetic energy operator in Eq. (24) and the chain rule (21) to go from Eq. (24) to (25). The final result (27) follows from Eq. (26) because $q_{t}$ is the expectation value of the position operator $\hat{q}$ in the state $\psi(t)$. Despite the energy conservation, the effective Hamiltonian determined by the effective potential of Eq. (23) is still time-dependent—the energy is conserved only because the Hamiltonian is nonlinear (i.e., it depends on the state $\psi$) and its change applied to $\psi$ happens to be “orthogonal” to the state $\psi$ [Eqs. (24)–(27)]. Therefore, the conservation of energy does not guarantee non-negative intensities in the spectrum. Yet, the hope is that the negative spectral features in the single-Hessian approach will be less pronounced than in the standard thawed Gaussian approximation.

The single-Hessian thawed Gaussian approximation may seem to be only a special, constant-Hessian case of one of several approximations used for accelerating semiclassical calculations based on the Herman-Kluk propagator.Herman and Kluk (1984) In Appendix A, we, therefore, compare the proposed method with the single-Hessian versions of the Herman–Kluk, Johnson’s,Gelabert et al. (2000) frozen Gaussian,Heller (1981b) adiabatic Herman–Kluk,Guallar, Batista, and Miller (1999) and prefactor-freeTatchen and Pollak (2009) approximations and show that the equivalence holds only for some of these methods and, moreover, only if the thawed Gaussian becomes “frozen,” which requires a specific choice of the reference Hessian.

II.6 Reference Hessians

Both global harmonic models and single-Hessian thawed Gaussian approximation depend on the choice of the reference Hessian. Two well-known special choices are the adiabatic Hessian—Hessian of the final electronic potential energy surface evaluated at its minimum ($V_{\text{ref}}=V_{\text{final}}$, $q_{\text{ref}}=q_{\text{eq,\ final}}$), and the vertical Hessian—Hessian of the final electronic surface evaluated at the Franck–Condon point, i.e., the minimum of the initial electronic surface ($V_{\text{ref}}=V_{\text{final}}$, $q_{\text{ref}}=q_{\text{eq,\ init}}$); see Fig. 2.Avila Ferrer and Santoro (2012); Egidi et al. (2014) We refer to the combinations of these two Hessian choices with the global harmonic approach as the adiabatic harmonic and vertical harmonic methods.Wehrle, Šulc, and Vaníček (2014); Patoz, Begušić, and Vaníček (2018); Begušić et al. (2018) In the literature, these global harmonic models are sometimes referred to as the adiabatic and vertical Hessian;Avila Ferrer and Santoro (2012); Egidi et al. (2014); Baiardi, Bloino, and Barone (2013) here, we use these names exclusively for the Hessians themselves to avoid the confusion between the single-Hessian thawed Gaussian propagation and global harmonic methods. The combinations of the single-Hessian approach with the different reference Hessians will be referred to as the adiabatic single-Hessian and vertical single-Hessian methods.

Finally, one can avoid computing any Hessian of the final electronic surface by using as reference the initial-state Hessian—Hessian of the initial electronic surface at its minimum ($V_{\text{ref}}=V_{\text{init}}$, $q_{\text{ref}}=q_{\text{eq,\ init}}$, see Fig. 2), which is commonly needed already for constructing the initial wavepacket. In the context of global harmonic methods, there are two natural possibilities of constructing a final-state harmonic potential using the initial-state Hessian: one can either compute the potential energy and gradient of the final-state potential energy surface at the initial geometry, which results in the vertical gradient model, or optimize the geometry in the final electronic state, which gives the adiabatic shift model.Avila Ferrer and Santoro (2012); Cerezo et al. (2013); Egidi et al. (2014) Both the vertical gradient and adiabatic shift models are examples of displaced harmonic systems, and thus ignore mode distortion and mixing (the Duschinsky effect) between the two electronic states. In the results section, we discuss only the adiabatic shift model and, for consistency with the other methods discussed in this work, refer to it as the initial harmonic model.

II.7 $PQ$-method and the Hamiltonian structure

The Riccati equation (11) can be solved with the “$PQ$ method”, i.e., by introducing auxiliary complex $D\times D$ matrices $Q_{t}$ and $P_{t}$ such thatHeller (1976a)

[TABLE]

Inverting Eq. (29) and inserting Eq. (28) into Eq. (11) yields the differential equations

[TABLE]

which can be recognized as Hamilton’s equations of motion

[TABLE]

of a “semiclassical” time-dependent Hamiltonian

[TABLE]

where $q_{t}$ plays a role of an external, time-dependent parameter. Above, ∗ denotes a complex conjugate and † the Hermitian transpose, i.e., a complex conjugate and transpose of a matrix. Hamilton’s equations (30)-(31) solve Eq. (11) for $\dot{A}_{t}$ for any choice of $Q_{0}$ and $P_{0}$ that satisfy Eq. (28) at time zero.

In the single-Hessian approximation, the time-dependent Hessian is replaced with the reference Hessian, and the semiclassical Hamiltonian (32) becomes independent of $q_{t}$ and, therefore, also independent of time. As a result, the quantum propagation using single-Hessian thawed Gaussian approximation for $H(q,p)$ can be mapped to exact classical propagation with a separable Hamiltonian

[TABLE]

Because of separability, both $H(q_{t},p_{t})$ and $H_{\text{sc}}(Q_{t},P_{t})$ are independent of time. In Appendix B, we show that the energy $E(t)$ of the wavepacket (7) is equal to $H_{\text{map}}(q_{t},p_{t},Q_{t},P_{t})$ for a specific choice of $Q_{0}$ and $P_{0}$ (up to a factor equal to Hagedorn parametrizationHagedorn (1980, 1998); Lubich (2008)), providing an independent proof of energy conservation by the single-Hessian thawed Gaussian approximation. Neither energy conservation nor mapping to a classical Hamiltonian system holds for the original thawed Gaussian approximation due to the dependence of the Hessian on $q_{t}$; in that case, Hamilton’s equation for $p_{t}$ derived from $H_{\text{map}}$ has an additional term compared to Eq. (10). Yet, a similar mapping, yielding a nonseparable Hamiltonian, does existKramer and Saraceno (1981); Arickx et al. (1986); Faou and Lubich (2006) if one applies the time-dependent variational principleHeller (1976b); Coalson and Karplus (1990) instead of the local harmonic approximation (8) to the quantum propagation of the Gaussian wavepacket (7).

III Computational details

III.1 Morse potential

To investigate the single-Hessian thawed Gaussian approximation in systems of varying anharmonicity, we constructed a series of Morse potentials,

[TABLE]

with different values of the dissociation energy $D_{e}$ and anharmonicity parameter $a$. In Eq. (34), $V_{\text{eq}}$ is the potential at the equilibrium position $q_{\text{eq}}$. We chose to work in atomic units ($\hbar=1$) and mass-scaled coordinates. The initial wavepacket was a real Gaussian with zero position and momentum, and with a width matrix $A_{0}=\omega_{0}/2$ corresponding to the ground vibrational state of a harmonic oscillator with frequency \omega_{0}=0.00456\ \text{a.u.}=1000\cm*-1*. The Morse parameters were $V_{\text{eq}}=0.1$ and $q_{\text{eq}}=\sqrt{2/\omega_{0}}=20.95\ \text{a.u.}$. We also fixed the global harmonic potential fitted to the Morse potentials at the equilibrium position $q_{\text{eq}}$; its frequency,

[TABLE]

was set to 0.0041\ \text{a.u.}=900\cm*-1*. Anharmonicity of the potential was controlled through the dimensionless constant

[TABLE]

Then, the $D_{e}$ and $a$ parameters were uniquely defined as

[TABLE]

The transition dipole moment was set to $1$. The wavepacket was always propagated for $4000$ steps of 8\a.u.\ \approx 0.194\fs. Spectra evaluated with the thawed-Gaussian, global harmonic, and single-Hessian approaches discussed in Section II.6 were compared with the exact quantum dynamics calculations, obtained with the second-order split-operator method. The position grid for the exact quantum dynamics consisted of $16384$ points between $-200$ and $200$ atomic units. To avoid artifacts of the finite-time calculation, all correlation functions were multiplied by a Gaussian damping function corresponding to the Gaussian broadening with half-width at half-maximum of 115\cm*-1*. Spectra were then computed from Eq. (4) and scaled according to the maximum intensity.

III.2 On-the-fly ab initio calculations

The on-the-fly ab initio implementation of the thawed Gaussian approximation has been detailed in Refs. Wehrle, Šulc, and Vaníček, 2014; Wehrle, Oberli, and Vaníček, 2015; Patoz, Begušić, and Vaníček, 2018; Begušić et al., 2018. Briefly, the method evaluates the required potential information along the trajectory from an ab initio electronic structure program. Our in-house code performs the dynamics, transformation between Cartesian and normal-mode coordinates, and interpolation of the Hessians if they are not computed at each step (see Refs. Wehrle, Šulc, and Vaníček, 2014; Patoz, Begušić, and Vaníček, 2018).

For ammonia, the ab initio calculations were performed using the complete active-space second-order perturbation theory, CASPT2(8/8), in combination with the aug-cc-pVTZ basis set, as implemented in the Molpro2012.1 package.Werner et al. (2012a, b) For the quinquethiophene, the ground-state potential data were evaluated using the density functional theory, while the time-dependent density functional theory was used for geometry optimization and Hessian calculations in the first excited electronic state; the functional was B3LYP and the basis set 6-31+G(d,p), as implemented in Gaussian09.Frisch et al. All trajectories were propagated using a time step of 8\a.u. for $1000$ steps in ammonia and for $997$ steps in quinquethiophene. In ammonia the Hessian was computed at each step, whereas in quinquethiophene the Hessian was evaluated only every four steps and interpolated in between; such an interpolation was previously validated in Ref. Wehrle, Šulc, and Vaníček, 2014. Before computing the spectra, the correlation functions were multiplied with a Gaussian damping function corresponding to the spectral Gaussian broadening with half-width at half-maximum of 200\cm*-1*. Further computational details about the ammonia absorption spectrum can be found in Ref. Wehrle, Oberli, and Vaníček, 2015 and, about the quinquethiophene emission spectrum, in Ref. Wehrle, Šulc, and Vaníček, 2014.

IV Results and discussion

IV.1 Morse potential

Figure 3 compares the exact spectra of two Morse potentials of different degrees of anharmonicity with those evaluated using the standard thawed Gaussian approximation, its adiabatic single-Hessian version, and the adiabatic harmonic method. In the weakly anharmonic potential (Fig. 3, left), all methods perform well, with only the global harmonic spectrum deviating slightly from the exact solution. In contrast, in the more anharmonic Morse potential, the adiabatic harmonic model recovers only the first few peaks. Interestingly, the single-Hessian version seems to be more accurate than the standard thawed Gaussian approximation in describing peak intensities.

To quantify the accuracy of the approaches discussed in Section II.6, we introduce the spectral contrast angle $\theta$ between a reference ($\sigma_{\text{ref}}$) and approximate ($\sigma$) spectra, conveniently defined through its cosine

[TABLE]

with the inner product $\sigma_{1}\cdot\sigma_{2}=\int d\omega\sigma_{1}(\omega)\sigma_{2}(\omega)$ of two spectra and norm of a spectrum $\lVert\sigma\rVert=\sqrt{\sigma\cdot\sigma}$. Spectra evaluated with the exact quantum dynamics are used as reference. In ab initio calculations, the errors in the absolute frequency shift of the spectrum originate mostly from the limited accuracy of the electronic structure methods. Therefore, even in the Morse potential, we first maximize the overlap with the reference by shifting the computed spectra in frequency and then evaluate the spectral contrast angle. The maximum overlap is found by scanning through all possible shifts, with the increment determined by the numerical resolution of the spectrum.

As shown in Fig. 4, the accuracies of all presented methods decrease with increasing anharmonicity of the potential. However, the methods based on the thawed Gaussian approximation clearly perform better than the global harmonic approaches. Moreover, the single-Hessian results are nearly the same for all three choices of the Hessian, which is not the case for the global harmonic approximations. The errors in the spectra of more anharmonic potentials (see Fig. 3) are reflected mainly in incorrect peak spacings, which are almost exclusively determined by the classical trajectory guiding the thawed Gaussian wavepacket—therefore, in the single-Hessian thawed Gaussian approximation the choice of the Hessian affects the result only weakly.

Negative intensities in the spectra computed with the thawed Gaussian approximation further increase the errors measured by the spectral contrast angle. Such features are nearly eliminated in the single-Hessian version of the thawed Gaussian approximation, which conserves energy exactly (see Fig. 5); however, negative intensities still arise even in the single-Hessian method due to the time dependence of the effective single-Hessian potential (23).

IV.2 Absorption spectrum of ammonia

Ammonia is a prototypical example of a floppy system, i.e., a system exhibiting large-amplitude motion. Electronic excitation to the first excited state is accompanied by a significant displacement of the umbrella inversion mode, allowing the generated wavepacket to visit anharmonic regions of the excited-state potential energy surface. Due to the small size of the system, rich nuclear dynamics, and available experimental data, the absorption, emission, and photoelectron spectra of ammonia have served as benchmarks for different methods built specifically to treat the anharmonicity effects.Tang, Imre, and Tannor (1990); Tang, Abramson, and Imre (1991); Capobianco et al. (2012); Baiardi, Bloino, and Barone (2017) In particular, the on-the-fly ab initio thawed Gaussian approximation showed significant improvement over the global harmonic models.Wehrle, Oberli, and Vaníček (2015)

Figure 6 compares the global harmonic and single-Hessian approaches with the on-the-fly ab initio thawed Gaussian approximationWehrle, Oberli, and Vaníček (2015) and with the experimental absorption spectrum of ammonia.Chen et al. (1999) All single-Hessian methods recover both the peak positions and intensities of the standard thawed Gaussian approximation. In contrast, all global harmonic models yield different and rather inaccurate spectra. Most interesting are the adiabatic single-Hessian thawed Gaussian approximation and adiabatic global harmonic model: although both methods use only one (adiabatic) Hessian, the former performs better than any other presented method, including the standard thawed Gaussian approximation, whereas the latter performs the worst. These results indicate that the single-Hessian thawed Gaussian approximation cannot be discarded in advance based on the performance of global harmonic models; in fact, its accuracy is much closer to that of the thawed Gaussian approximation. Indeed, even the initial (ground-state) single Hessian approach reproduces almost perfectly the result of the standard on-the-fly ab initio thawed Gaussian approximation.

IV.3 Emission spectrum of

quinquethiophene

Due to their potential in molecular electronics, polythiophenes and their derivatives have been studied extensively. Oligothiophenes have also served as a model system for studying the dependence of optical properties on the system size. They present a challenge for computing vibrationally resolved electronic spectra due to the torsional degrees of freedom, which cannot be treated with global harmonic models. Wehrle et al.Wehrle, Šulc, and Vaníček (2014) showed that the on-the-fly ab initio thawed Gaussian approximation performs well despite the double-well character of the potential along the torsional modes connecting the planar and twisted structures.

In Fig. 7, we compare the experimentalBecker et al. (1996) emission spectrum of quinquethiophene, an oligomer composed of five thiophene units, and corresponding spectra computed with various approximations discussed in Section II. The single-Hessian approaches using the initial (excited-state) and vertical Hessians produce almost the same spectra as the standard thawed Gaussian approximationWehrle, Šulc, and Vaníček (2014) (shown in Fig. 7, top). However, this is not the case for the adiabatic single-Hessian method, which yields a broad spectrum due to the incorrect description of the torsional degrees of freedom. As discussed in Ref. Wehrle, Šulc, and Vaníček, 2014, the initial wavepacket is placed at the top of a potential barrier along the torsional modes, which results in a constant but slow wavepacket spreading. The adiabatic Hessian has all frequencies positive and is therefore qualitatively inappropriate. Interestingly, the initial single-Hessian approach, which propagates a frozen Gaussian, results in a rather accurate spectrum, implying that the errors of using the adiabatic Hessian arise due to the incorrect width of the Gaussian wavepacket.

In contrast, the failure of the adiabatic global harmonic model (Fig. 7, top right) is not related to the Hessian, but rather to the large displacement of the ground-state potential minimum from the initial geometry. The computed emission spectrum is nearly featureless because the wavepacket quickly drifts away from the initial planar geometry and does not return during the short dynamics considered for spectra simulations. This explanation is supported by the equally featureless spectrum of the adiabatic shift model, i.e., the initial harmonic model (see Fig. 7, bottom right), which has the same displacement of the adiabatic global harmonic model but uses the initial (excited-state) Hessian.

V Conclusion

In conclusion, we have presented and validated an efficient method for evaluating low-resolution vibronic spectra of polyatomic molecules. The proposed single-Hessian thawed Gaussian approximation, whose computational cost lies between those of the global harmonic and thawed Gaussian approximations, performs surprisingly well, in some cases even better than the more computationally demanding thawed Gaussian approximation. Moreover, unlike the standard thawed Gaussian approximation, the single-Hessian approach conserves total energy exactly. We have shown that despite the conservation of energy, the computed spectra may still contain negative intensities due to the time dependence of the effective Hamiltonian. Yet, the negative spectral features are significantly smaller compared with the standard thawed Gaussian approximation. In contrast to the spectra evaluated using the global harmonic approaches, those computed with the single-Hessian thawed Gaussian approximation depend only weakly on the reference Hessian. Therefore, the single-Hessian approach offers a considerable and systematic improvement over the commonly used global harmonic models at the cost of a single ab initio classical trajectory.

Acknowledgements.

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

Appendix A Single-Hessian approximations of the Herman–Kluk

prefactor

Within the Herman–KlukHerman and Kluk (1984); Herman (1986); Kay (1994a, b) semiclassical initial value representation,Miller (1970, 2001); Kay (2005) the quantum evolution operator is approximated as

[TABLE]

where $D$ is the number of degrees of freedom, $S_{t}:=\int_{0}^{t}L_{t^{\prime}}dt^{\prime}$ the classical action,

[TABLE]

the Herman–Kluk prefactor, $M_{ab}=\partial a_{t}/\partial b_{0}$ components of the stability matrix, $|q_{t}p_{t}\rangle$ the coherent state whose wavefunction in position representation is

[TABLE]

$\Gamma$ denotes a pure imaginary symmetric coherent state width matrix (i.e., $\Gamma^{\ast}=-\Gamma$ and $\Gamma^{T}=\Gamma$) , and $q_{t}$ and $p_{t}$ evolve classically according to Eqs. (9)–(10).

Reversing the main idea of the $PQ$ methodHeller (1976b) mentioned in Section II.7, in the log-derivative formulation,Gelabert et al. (2000); Di Liberto and Ceotto (2016) the Herman–Kluk prefactor is expressed in terms of an auxiliary matrix

[TABLE]

where $Q_{t}=M_{qq}\cdot Q_{0}+M_{qp}\cdot P_{0}=M_{qq}-M_{qp}\cdot\Gamma$, as

[TABLE]

Matrix $Q_{t}$ defined here is equivalent to that of Eq. (28) for a specific choice of initial conditions: $Q_{0}=I_{D}$ and $P_{0}=-\Gamma$ (see Appendix B). Matrix $\alpha_{t}$, whose initial value is $\alpha_{0}=-\Gamma$, obeys the same equation of motion as the matrix $A_{t}$ of the thawed Gaussian approximation [Eq. (11)]; the connection between $\alpha_{t}$ and $A_{t}$ was discussed, e.g., in Ref. Gelabert et al., 2000.

To compare different approximations to the prefactor $R_{t}$ with the single-trajectory single-Hessian thawed Gaussian approximation, we consider only a single trajectory in Eq. (40) and approximate the propagated wavepacket as

[TABLE]

Then, $\Gamma=A_{0}$ of the initial wavepacket $\psi_{i}$ and the wavepacket at time $t$ is a Gaussian (7) with parameters $A_{t}=A_{0}$ and $\gamma_{t}$ given by

[TABLE]

In what follows, we apply the single-Hessian potential [Eq. (23)] to the Herman–Kluk prefactor and its approximations. For a constant Hessian, assuming for simplicity that $\Gamma$, $m$, and $V_{\text{ref}}^{\prime\prime}(q_{\text{ref}})$ commute (which is valid, e.g., if $D=1$, or if all three matrices are diagonal, or if spherical Gaussians and mass-scaled coordinates are used, i.e., $\Gamma\propto I_{D}$ and $m\propto I_{D}$), the Herman–Kluk prefactor simplifies toGelabert et al. (2000)

[TABLE]

Matrix $\tilde{\alpha}_{t}$ evolves as $A_{t}$ of the single-Hessian thawed Gaussian approximation, but with a modified initial condition

[TABLE]

where $A_{\text{ref}}$ corresponds to the coherent state of a harmonic potential with force constant matrix $k=V_{\text{ref}}^{\prime\prime}(q_{\text{ref}})$, i.e.,

[TABLE]

Equation (47) coincides with the slowly varying Hessian approximation of Gelabert et al.Gelabert et al. (2000); however, their approximation formally assumes a time-dependent Hessian for the evolution of $\tilde{\alpha}_{t}$, whereas here, Eq. (47) is an exact expression for the prefactor in the approximate potential (23). Because matrix $\tilde{\alpha}_{t}$ is, in general, complex at $t>0$, the norm of the “single-Hessian Herman–Kluk” wavepacket is not conserved. This is remedied easily by taking only the imaginary part of $\tilde{\alpha}_{t}$ in Eq. (47), or, equivalently, by renormalizing the wavepacket at each step.

In Johnson’s multichannel Wentzel–Kramers–Brillouin approximation,Gelabert et al. (2000); Issack and Roy (2005, 2007a, 2007b) one assumes that the matrix $\alpha_{t}$ varies slowly, i.e., $\dot{\alpha}_{t}\approx 0$, which yields

[TABLE]

where $\omega_{j}(t)$ are time-dependent frequencies obtained from the Hessian evaluated at $q_{t}$. This method involves a time-dependent Hessian and is, therefore, closer to the original thawed Gaussian approximation than to its single-Hessian version. However, if Johnson’s approximation is combined with the single-Hessian potential (23), the time-dependent frequencies $\omega_{j}(t)$ are replaced by the reference frequencies $\omega_{\text{ref},j}$ obtained from the reference Hessian $V_{\text{ref}}^{\prime\prime}(q_{\text{ref}})$ and the integral in Eq. (51) is trivial.

The adiabatic approximationGuallar, Batista, and Miller (1999, 2000) of the Herman–Kluk prefactor assumes an instantaneously diagonal Hessian at each time step, i.e., it neglects the offdiagonal entries of the full Hessian matrix. Within the single-Hessian approximation, the resulting expression for $R_{t}$ is the same as for the single-Hessian Herman–Kluk [Eq. (47)] except for a modified (diagonal) Hessian.

Finally, the crudest approximation is to replace the prefactor by unity, which is known as the prefactor-free approach;Tatchen and Pollak (2009) then, $\gamma_{t}=S_{t}$ and no Hessian computation is needed.

Equations of motion for parameters $A_{t}$ and $\gamma_{t}$ in the single-Hessian thawed Gaussian, Herman–Kluk, Johnson’s, adiabatic Herman–Kluk, and prefactor-free approximations are summarized in Table 1, where we also present analogous expressions for Heller’s frozen Gaussian approximation.Heller (1981b) Single-Hessian thawed Gaussian wavepacket has a time-dependent width, whereas the other approximations propagate a coherent state with only a modified phase factor. Special case is the initial single-Hessian approach, which uses the initial Hessian for the single-Hessian thawed Gaussian propagation. Then, $\dot{A}_{t}=0$ holds even for the thawed Gaussian wavepacket and $\gamma_{t}$ is the same for the thawed Gaussian, Herman–Kluk, Johnson’s, and frozen Gaussian approximations. Let us emphasize that in the multiple-trajectory implementations of the single-Hessian Herman–Kluk, Johnson’s, and frozen Gaussian methods, $\Gamma$ is a free parameter; for $\Gamma=A_{\text{ref}}$, the three approximations are equivalent. In contrast, in the single-trajectory thawed Gaussian approximation, because the initial width parameter $A_{0}$ is fixed by the wavepacket $\psi_{i}$, using a constant Hessian does not imply a time-independent matrix $A_{t}$. Therefore, the single-Hessian method is, despite similarities, fundamentally different from other approaches.

Various single-Hessian approaches are compared numerically in Table 2. The results confirm that the single-Hessian thawed Gaussian approximation is not identical to the single-trajectory Herman-Kluk propagator or any of its several simplified versions. That the differences between the methods are only small may be attributed to a weak distortion of the model system—greater difference between the ground- and excited-state Hessians would lead to greater deformations of the wavepacket, which cannot be described by a single coherent state [Eq. (45)]. The shifted spectra obtained with single-Hessian Johnson’s, frozen Gaussian, and prefactor-free approximations are the same because the methods differ only by a factor $\exp(it\Delta)$, where $\Delta$ is a real constant depending on the methods that are compared (see Table 1). Finally, all methods except for the prefactor-free approximation yield exactly the same result if the initial-state Hessian is used as a reference, in agreement with the theoretical justification given above.

The single-Hessian approximations of the coherent-state methods are not necessarily useful in practice and are presented here only for comparison with the single-Hessian thawed Gaussian approximation. Indeed, in the usual multi-trajectory setup, the single-Hessian Herman–Kluk approach, which is equivalent to the harmonic approximation mentioned briefly in Ref. Di Liberto and Ceotto, 2016, would already be a feasible computational method and no further approximations of the prefactor would be needed. Otherwise, approaches based on the Herman–KlukIanconescu, Tatchen, and Pollak (2013); Bonfanti et al. (2018) and Johnson’sIssack and Roy (2005, 2007a, 2007b) approximations have been validated on difficult systems, where accurate calculations require the evaluation of Hessians along each trajectory.

Appendix B Energy of the Gaussian wavepacket and the

mapping Hamiltonian

B.1 Useful relations

Auxiliary matrices $Q_{t}$ and $P_{t}$, defined by Eqs. (28) and (29), satisfy the relationsLee and Heller (1982); Faou, Gradinaru, and Lubich (2009)

[TABLE]

The former is obtained from $P_{t}^{T}\cdot Q_{t}=(Q_{t}^{T}\cdot P_{t})^{T}$ using Eq. (28) for $P_{t}$, the latter by showing that the time derivative of the left hand side is zero and by confirming the relation at time zero—by realizing that

[TABLE]

A remarkable relationHagedorn (1998); Faou, Gradinaru, and Lubich (2009)

[TABLE]

can be deduced from Eqs. (52) and (53):

[TABLE]

Equation (57) follows from Eq. (56) because $Q_{t}^{\dagger}\cdot P_{t}^{\ast}=(Q_{t}^{T}\cdot P_{t})^{\ast}=(P_{t}^{T}\cdot Q_{t})^{\ast}=P_{t}^{\dagger}\cdot Q_{t}^{\ast}$, where we used Eq. (52), and Eq. (58) follows from (57) by applying Eq. (53).

B.2 Energy of the thawed Gaussian wavepacket

The total energy of the thawed Gaussian wavepacket, computed as the expectation value $E=\langle\psi(t)|\hat{H}_{\text{eff}}(t)|\psi(t)\rangle$, can be split asWehrle, Oberli, and Vaníček (2015)

[TABLE]

into the “classical” energy of the central trajectory,

[TABLE]

and “semiclassical” energy

[TABLE]

The first factor inside the trace can be rewritten as

[TABLE]

Equation (62) holds because $A_{t}$ is symmetric, in Eq. (63) we used expression (28) for $A_{t}$, and in Eq. (64) we introduced a matrix-valued function

[TABLE]

The last step (65) follows because both $E_{\text{sc}}$ and $\operatorname{Im}A_{t}$ in Eq. (61) are real. As for the second factor inside the trace in Eq. (61), relations (28) and (54) imply that

[TABLE]

Substitution of expressions (67) and (64) for the two factors into the relation (61) for the semiclassical energy gives

[TABLE]

The choice of $Q_{0}$ is not determined by the definitions (28) and (29) of $Q_{t}$ and $P_{t}$. A common choice is $Q_{0}=I_{D}$ (a $D$-dimensional identity matrix) and $P_{0}=A_{0}$, which yields

[TABLE]

However, one can remove all constant factors from Eq. (68) by setting $Q_{0}=(2\operatorname{Im}A_{0}/\hbar)^{-1/2}\cdot U$, with an arbitrary unitary matrix $U$, to obtain

[TABLE]

where $H_{\text{sc}}$ is the semiclassical Hamiltonian (32) from Section II.7. Note that with this choice of $Q_{0}$, the right-hand side of the generalized commutation relation (53) becomes $i\hbar I_{D}$, in direct analogy with $[\hat{q},\hat{p}]=\hat{q}\otimes\hat{p}^{T}-\hat{p}\otimes\hat{q}^{T}=i\hbar I_{D}$, but differs slightly from Hagedorn’s conventionHagedorn (1980, 1998) of $2iI_{D}$, which would also fail to eliminate the factor $\hbar/2$ in the energy (69). If the exact potential is replaced with the single-Hessian potential $V_{\text{SH}}$ [Eq. (23)], the matrix function $\mathcal{H}_{\text{sc}}$ from Eq. (66) becomes independent of $q_{t}$, and so does $H_{\text{sc}}$. As discussed in Section II.7, in this setting $H_{\text{sc}}(Q_{t},P_{t})$ is a constant of motion, and, so is the semiclassical energy, since, according to Eq. (70), it is equal to the semiclassical Hamiltonian $H_{\text{sc}}$. Finally, in agreement with the derivation presented in Sec. II.5, the total energy is conserved because it is equal to the mapping Hamiltonian $H_{\text{map}}(q_{t},p_{t},Q_{t},P_{t})$ of Eq. (33).