On-the-fly ab initio semiclassical evaluation of vibronic spectra at finite temperature

Vibrationally resolved electronic spectra have, for a very long time, been used to learn more about electronic and vibrational states of molecules, their potential energy surfaces, and light-induced dynamics of nuclei.^1–4 The computational methods for simulating such spectra are, therefore, an essential tool in physical chemistry.

The most widespread is the global harmonic method,^5–7 which employs the harmonic approximation for both ground- and excited-state potential energy surfaces. Within the framework of the global harmonic approximation, one can easily account for non-Condon and finite-temperature effects.^8–12 This approximation, however, neglects the effects of anharmonicity, which can significantly alter molecular spectra. Other quantum^13–16 and semiclassical^4,15,17–20 methods do include anharmonicity effects on spectra, but at a substantial computational cost. Recently, we have been investigating the thawed Gaussian approximation (TGA),^21–24 an efficient semiclassical method that accounts partially for anharmonicity and requires no initial knowledge of the potential energy surface. The method has been extended to include non-Condon effects, namely, to account for the Herzberg–Teller contribution to the transition dipole moment.^25–27 Unfortunately, as a wavepacket propagation method, it has been limited to computing spectra in the zero-temperature limit, where only the ground vibrational state is populated initially.

To account for non-zero temperature, one typically employs the density matrix formalism, where a number of numerically exact^28–30 and approximate^31–37 approaches exist. Otherwise, typical wavefunction-based methods can be used in combination with statistical sampling of initial conditions.^38–44 Thermo-field dynamics^45,46 offers an alternative way to make wavefunction-based methods applicable at finite temperature: the problem, which seemingly requires the von Neumann equation for the density matrix, is mapped to a time-dependent Schrödinger equation with twice as many degrees of freedom. Recently, the thermo-field dynamics was employed in chemistry for solving the coupled electronic-vibrational dynamics,^47–50 electronic structure,⁵¹ and vibronic spectroscopy¹¹ problems. The application to vibronic spectroscopy, which is of central interest to this work, was, however, restricted to the global harmonic approximation.

Here, we combine the extended thawed Gaussian wavepacket propagation with the thermo-field dynamics in order to include both anharmonicity and finite-temperature effects. Due to the favorable scaling of the thawed Gaussian approximation with the system’s size, the new method adds nearly no additional cost to the original, zero-temperature approach. To illustrate the accuracy achieved by going beyond both global harmonic and zero-temperature approximations, we test the method on a set of Morse potentials with different degrees of anharmonicity and at different temperatures. Finally, we apply it to evaluate the spectrum corresponding to the symmetry-forbidden electronic transition S₁ ← S₀ (⁠$Ã1B2u←X̃1A1g$⁠) of benzene and demonstrate that the simultaneous inclusion of Herzberg–Teller, anharmonicity, and finite-temperature effects is needed to reproduce the experimental spectrum.

Before turning to vibrationally resolved electronic spectra at finite temperature, let us briefly describe the original, zero-temperature approach based on the extended thawed Gaussian approximation.

The absorption spectrum at zero temperature can be computed as the Fourier transform,^33,52,53

of the correlation function,

where |1, g⟩ is the ground (“g”) vibrational state of the ground (“1”) electronic state, ℏω_1,g is its energy, Ĥ₂ is the nuclear Hamiltonian corresponding to the excited (subscript “2”) electronic state, and $μ^$ is the transition dipole moment $μ^21=μ→^21⋅ϵ→$ projected on the polarization $ϵ→$ of the external electric field. In other words, to compute the spectrum, one has to evolve the nuclear wavefunction $|ϕ0=μ^|1,g$ on the excited-state surface, which is, in general, a challenging task that scales exponentially with the number of atoms.

Different exact quantum^42,54–57 and semiclassical^19,58–74 methods were developed for solving the problem of wavepacket propagation. Sometimes, the region of the excited-state potential energy surface explored by the evolved wavepacket is fairly harmonic, meaning that it can be well approximated by a second-order Taylor expansion in nuclear coordinates about a reference geometry; we call this the global harmonic approximation. Then, the correlation function (2) can be obtained analytically.^8–10,75,76 Moreover, in this case, the excited-state surface is easily constructed from a single Hessian calculation. The prevalence of the global harmonic method in the vibronic spectroscopy literature^{7,8,10,75–79} testifies, on the one hand, to its applicability in a wide range of molecules and, on the other hand, to the absence of accessible alternatives that can account for anharmonicity effects. Often, due to the reduced resolution of electronic spectra, even such a crude approximation, which would nowadays be almost unacceptable for the simulation of vibrational (infrared) spectra, is considered appropriate.

To account for the anharmonicity effects on the spectrum at least approximately, we recommend using the simple and efficient semiclassical thawed Gaussian approximation.²¹ In contrast to many other exact or approximate quantum dynamics methods, this method is computationally feasible even for rather large molecules and can be employed in a “black-box” fashion, i.e., it requires little human input. In particular, the thawed Gaussian approximation requires only local potential energy information along the classical trajectory (as described below) and, therefore, supports an on-the-fly implementation where the potential energy is provided by an ab initio electronic structure calculation.

Within the thawed Gaussian approximation,²¹ the wavepacket is assumed to be a complex Gaussian function,

parameterized by the time-dependent D-dimensional real vectors q_t and p_t, D × D complex symmetric matrix A_t, and complex number γ_t; D is the number of degrees of freedom. The time dependence of the matrix A_t implies that the width of the thawed Gaussian wavepacket changes with time, as opposed to the frozen Gaussian ansatz where the width remains constant. Wavepacket (3) solves exactly the Schrödinger equation,

where $T(p)=12pT⋅m−1⋅p$ is the kinetic energy and

is the local harmonic approximation of the true potential energy V(q) about q_t, if the time-dependent parameters of ψ_t satisfy the following equations of motion:²¹ 

In these equations, V′(q_t) and V″(q_t) denote, respectively, the gradient and Hessian of the potential energy evaluated at q_t, m is the symmetric mass matrix, and L_t = T(p_t) − V(q_t) is the Lagrangian. Note that due to the local harmonic approximation, Eq. (4) is a nonlinear Schrödinger equation because the potential V_LHA(q) depends on the wavefunction through the parameter q_t [Eq. (5)], i.e., V_LHA(q) ≡ V_LHA(q; q_t) ≡ V_LHA(q; ψ_t).

To construct the initial wavepacket, the ground-state potential energy surface V₁(q) is assumed to be harmonic in the vicinity of its minimum q_eq, i.e., the ground-state Hamiltonian is approximated as

where $K=V1′′(qeq)$ is the symmetric force-constant matrix and ∂_q = ∂/∂q. In position representation, the lowest eigenstate ψ₀(q) of the Hamiltonian (10) is a Gaussian (3) with parameters

where $Ω=m−1/2⋅K⋅m−1/2$⁠. This initial wavefunction ψ₀(q) is then evolved by solving differential equations (6)–(9) with V = V₂ (the excited-state potential energy).

The thawed Gaussian wavepacket (3) is not suited to treat non-Condon effects, i.e., the effects due to the dependence of the transition dipole moment μ(q) on nuclear coordinates q. Within the Herzberg–Teller approximation—the simplest extension of the Condon approximation—the transition dipole moment is assumed to be a linear function,⁸⁰ 

where μ′(q_eq) is the gradient of μ with respect to nuclear coordinates at the equilibrium geometry. Then, ϕ₀(q) = μ(q)ψ₀(q) is no longer a Gaussian wavepacket. Fortunately, the extended thawed Gaussian ansatz,^25,81

which is a special case of Hagedorn’s “Gaussian times a polynomial” wavepacket,^82–84 solves the same Schrödinger equation [Eq. (4)] as ψ_t(q), provided that the Gaussian parameters evolve, as before, according to Eqs. (6)–(9) and, in addition,

Hence, with the extended thawed Gaussian approximation, one can include the Herzberg–Teller contribution at nearly no additional computational cost.

At non-zero temperature, the dipole–dipole correlation function needed in vibronic spectroscopy is

where $ρ^=e−βĤ1/Tr(e−βĤ1)$ is the vibrational density operator and β = 1/k_BT. Note that in Eq. (19), we assumed that only the ground electronic state is populated in the thermal equilibrium, which is usually justified by the large energy gap between the ground and first excited electronic states. Because the time evolution in Eq. (19) involves two different Hamiltonians, an obvious classical analog of Eq. (19) is missing, which, in turn, hinders the development of classical-like or semiclassical approximations for C(t). Here, we demonstrate that by transforming the problem to the one of wavepacket dynamics in an augmented space, one can easily make use of the existing semiclassical methods for solving the time-dependent Schrödinger equation.

The correlation function can be re-written as

where

$q¯=(q,q′)T$is a 2D-dimensional coordinate vector, and

is a Hamiltonian in $q¯$ coordinates. In Eq. (20), we used the relation $[ρ^,Ĥ1]=0$ and the cyclic property of the trace; in Eq. (21), we introduced the position representation in q and q′ coordinates; in going from (25) to (26), we used the fact that the two Hamiltonians H₁(q′) and H₂(q) commute because they act on different coordinates; finally, Eq. (22) follows from (21) because

and $(ρ^1/2)†=ρ^1/2$⁠.

Equation (22) has a remarkable interpretation—the dipole–dipole correlation function C(t) for a D-dimensional system at finite temperature T can be thought of as a wavepacket autocorrelation function of $ϕ¯t(q¯)$ evolved with the Hamiltonian $H¯(q¯)$ according to the Schrödinger equation,

which describes an effective 2D-dimensional system at zero temperature.

The approach described here is, despite the explicit use of the position representation, equivalent to the thermo-field dynamics, as presented in Ref. 11. Indeed, the final result does not depend on the representation,

where $|q¯=|q|q̃′$ is a general position state in the augmented direct-product Hilbert space and $|q̃$ denotes a position state in the “fictitious” (or “tilde”) Hilbert space. In  Appendix A, we derive Eq. (31) using the standard thermo-field dynamics notation and without invoking the position representation.

In principle, any known method for solving the time-dependent Schrödinger equation can be applied to obtain $ϕ¯t$⁠. However, the doubled number of coordinates adds a substantial, if not prohibitive, computational cost to the already large cost of zero-temperature calculations with exponentially scaling exact quantum methods. In the following, we therefore employ the extended thawed Gaussian approximation, which scales favorably with the number of degrees of freedom.

To solve Eq. (30) with the (extended) thawed Gaussian approximation, we must first identify $ϕ¯0$ and the local harmonic approximation to the potential energy $V¯(q¯)=V2(q)−V1(q′)$⁠.

If we assume, as in Sec. II A, that the ground-state surface V₁ is harmonic [Eq. (10)], a general off-diagonal matrix element $ρ1/2(q,q′)≡ρ1/2(q¯)$ of $ρ^1/2$ is a Gaussian (3) parameterized with 2D-dimensional vectors,

a 2D × 2D matrix,

composed of D × D submatrices

and a scalar

See  Appendix B for the derivation of Eqs. (32)–(35). The harmonic approximation for the ground-state potential energy surface is justified in the vicinity of its minimum and, therefore, for the construction of the equilibrium vibrational density matrix. In fact, even in fairly anharmonic systems, the Gaussian density matrix often serves as a good starting point for semiclassical approximations.^85–88 Next, we assume $μ^$ to be diagonal in position representation,

and employ the Herzberg–Teller approximation [Eq. (15)] to obtain

With these initial values, we propagate the time-dependent parameters $q¯t$⁠, $p¯t$⁠, Ā_t, and $γ¯t$ according to Eqs. (6)–(9) and $b¯t$ according to Eq. (18). The potential energy, its gradient, and its Hessian are given by

while the D × D mass matrix m is replaced by the 2D × 2D matrix,

where q_t and $qt′$ are D-dimensional vectors composed of the first and second halves of coordinates of $q¯t$⁠, i.e., $q¯t=(qt,qt′)T$⁠. Interestingly, the classical equations of motion [Eqs. (6) and (7)] for the parameters $q¯t$ and $p¯t$ are solved by propagating two independent trajectories in D spatial dimensions: the first trajectory evolves q_t and p_t with the excited-state Hamiltonian H₂, while the second trajectory evolves $qt′$ and $pt′$ with the negative of the ground-state Hamiltonian, −H₁, due to the negative signs of mass in Eq. (42) and gradient in Eq. (40). Because the second trajectory is at a fixed point, i.e., at the minimum of the ground-state potential energy V₁ with zero momentum, it shows no dynamics. As a result, one requires only a single excited-state classical trajectory, to evolve the first D coordinates of $q¯t$⁠, and Hessians of the excited-state potential energy surface along this trajectory, which is the same as in the original zero-temperature approach; no further potential energy evaluations are needed to account for the temperature effects.

Finally, let us note that an alternative approach to finite-temperature spectra with Gaussian wavepackets has been proposed in Ref. 12. There, the authors proposed a similar scheme to directly evolve the coherence $ρ^μ(t)=exp(−iĤ2t/ℏ)μ^ρ^exp(iĤ1t/ℏ)$ in position representation with the doubled number of degrees of freedom. Then, the correlation function is evaluated simply as $C(t)=Tr[μ^†ρ^μ(t)]$⁠. Their method, combined with the local harmonic approximation, is equivalent to ours and gives the same correlation function. In contrast to our approach, the method of Ref. 12 has so far been used to compute vibronic spectra only in systems described with globally harmonic potential energy surfaces, where it is equivalent to the global harmonic approximation for vibronic spectra, which is much simpler because analytical expressions for C(t) exist.^10,75 Our approach based on thermo-field dynamics has the advantage of reducing the transition dipole autocorrelation function to the simpler and well-known expression (31) for the wavepacket autocorrelation, thus making it very easy to implement the finite-temperature treatment of vibronic spectra into the standard zero-temperature wavefunction-based codes, which typically contain procedures for computing the wavepacket autocorrelation.

To test the accuracy of the proposed method, we construct a one-dimensional model system consisting of a ground-state harmonic potential and an excited-state Morse potential. The ground-state surface is assumed harmonic to exclude the error (or error cancellation) due to using an approximate initial vibrational state—this is rarely an issue with zero-temperature methods because the harmonic approximation typically holds in the vicinity of the potential minimum but could affect the results at higher temperatures. In the current model, the error of the results obtained with the thawed Gaussian approximation is only due to the anharmonicity of the excited-state potential energy surface.

A set of Morse potentials,

was constructed by fixing the equilibrium position q₂, minimum energy V₂(q₂), and frequency $ω2=V2′′(q2)/m$ at q₂ and by varying the anharmonicity parameter χ. We set the minimum of the ground-state harmonic potential to zero (q₁ = 0) and its frequency to ω₁ = 1, while the excited-state Morse parameters were q₂ = 1.5, ω₂ = 0.9, and V(q₂) = 10. Mass was set to m = 1. The level of anharmonicity was tuned by changing the parameter χ in the range between 0.01 and 0.02 in steps of 0.001.

The exact spectrum was computed by evaluating Franck–Condon factors by numerical integration, which is feasible for this one-dimensional model system since both harmonic and Morse vibrational eigenfunctions are known analytically. The adiabatic harmonic model,

which is constructed about the minimum of the potential energy surface, is the same for all constructed Morse potentials because it does not depend on χ. Since the (extended) thawed Gaussian approximation is exact for harmonic potentials, it was used to compute the adiabatic harmonic spectra. For both harmonic and thawed Gaussian dynamics calculations, the time step was 0.1 and the total simulation time was 1000, i.e., 10 000 steps in total were taken. Gaussian broadening with the half-width at half-maximum of 0.1 was applied to all spectra. Spectra were evaluated at scaled temperatures T_ω = 0, 0.5, and 1, where T_ω = k_BT/ℏω₁ = 1/βℏω₁ (e.g., for an average molecular vibration of ω = 1000 cm⁻¹, T_ω = 1 corresponds to the temperature T ≈ 1439 K). A constant transition dipole moment μ = 1 was used.

To compare reference (σ_ref) and approximate (σ) spectra, we used the spectral contrast angle θ, defined through its cosine as

where σ₁ · σ₂ = ∫dωσ₁(ω)σ₂(ω) is the inner product of two spectra and $∥σ∥=σ⋅σ$ is the associated norm. In all calculations, the reference was the exact spectrum, while the approximate spectra were computed with the adiabatic global harmonic and thawed Gaussian approximations.

The S₁ ← S₀ absorption spectrum of benzene was computed with adiabatic harmonic, vertical harmonic, and thawed Gaussian approximations. In short, the adiabatic harmonic model is, as described above, obtained by the second-order Taylor expansion of the excited-state potential energy surface about its minimum, while for the vertical harmonic model, the same expansion is performed about the ground-state minimum.

Density functional theory was used for the optimization and Hessian calculation of the ground electronic state, while its time-dependent version was employed for the excited-state optimization, energy, gradient, and Hessian calculations. We used the B3LYP functional with the ultrafine grid and 6-31+G(d,p) basis set, as implemented in the Gaussian09⁸⁹ package. For the thawed Gaussian propagation, we used a second-order symplectic integrator with a time step of 8 a.u. (≈0.2 fs) and 10 000 steps in total. The Hessian of the potential energy was evaluated every four steps and interpolated in between, as done previously in Ref. 25. The ground-state surface was assumed to be harmonic. The gradient of the electronic transition dipole moment was computed numerically by the second-order finite difference method with a step of 10⁻⁴ Å.²⁵ 

The computed correlation functions were multiplied by an exponential damping function e^−t/τ with τ = 18 000 a.u., resulting in a Lorentzian line shape with the half-width at half-maximum of ≈12.2 cm⁻¹. To facilitate comparison with the experimental spectrum, computed spectra were shifted and scaled to match the experimental spectrum of Ref. 90 at its highest peak (data taken from the MPI-Mainz UV/VIS Spectral Atlas^91,92).

Finally, let us emphasize that the finite-temperature treatment of spectra requires no additional electronic structure evaluations, i.e., the same ab initio data could be reused to compute the benzene spectrum at any given temperature. We evaluated the spectra at zero temperature and at the temperature of the experiment (T = 298 K).

Thawed Gaussian and global harmonic spectra were compared with the exact result (see Fig. 1). Already for the system with weak anharmonicity (left panels, χ = 0.01), the thawed Gaussian approximation provides a more accurate spectrum than the harmonic method. The difference is seen mainly in the intensities of the high-frequency peaks. Since the adiabatic harmonic model describes well the region around the potential minimum, it can recover the positions and intensities of peaks corresponding to transitions between vibrational states with small quantum numbers. In contrast, the harmonic approximation breaks down for vibrational states with more quanta, resulting in incorrect intensities of high-frequency transitions. The effect of anharmonicity on the peak positions becomes significant for χ = 0.02, and even the thawed Gaussian approximation is inadequate. Nevertheless, it is still more accurate than the adiabatic harmonic model, which has no dependence on χ (harmonic spectra are, clearly, the same for different χ at a given temperature). In typical molecular systems, the peaks are often left unresolved due to the short excited-state lifetime or inhomogeneous broadening. Then, the intensities play an important role in recovering the overall shape of the spectrum, whereas even an error of tens of reciprocal centimeters in peak positions can be tolerated.

In contrast to the global harmonic method, the thawed Gaussian approximation can result in non-physical negative spectral features, which are due to the nonlinear character of the Schrödinger equation (4). This is a well-known disadvantage of the method and was discussed in more detail elsewhere.^23,24 In the studied Morse system, a negative peak overlaps with the hot band around ω = 9, resulting in a poor description of this spectral region at higher temperatures. In a way, the gain in accuracy in the high-frequency part of the spectrum is accompanied by a loss in accuracy in the frequency region below the 0–0 transition.

To compare the global harmonic and thawed Gaussian methods quantitatively, we measure the error of an approximate spectrum with the spectral contrast angle between the approximate and exact spectra (see Fig. 2). The thawed Gaussian approximation gives more accurate spectra than the harmonic approximation for all anharmonicities and at all temperatures studied. However, an interesting trend is observed: the harmonic approximation becomes more accurate as the temperature increases, whereas the thawed Gaussian approximation keeps the same degree of accuracy at all temperatures. The main reason for such behavior is closely related to the discussion above. As the temperature increases, the intensity of hot bands below the 0–0 transition grows, and they become more relevant in measuring the error. Hence, the adiabatic harmonic method gains on accuracy, unlike the thawed Gaussian approximation, which always loses on accuracy in the low-frequency part of the spectrum. However, such behavior of the global harmonic method is not general; if the ground-state potential energy surface were anharmonic, high-temperature spectra would also reflect the effects neglected in the global harmonic models—those of ground-state anharmonicity on the initial density matrix.

The symmetry-forbidden S₁ ← S₀ transition in benzene is a well-known example of the Herzberg–Teller effect,^1,2 where the spectrum arises only due to the coordinate dependence of the transition dipole moment, which is zero by symmetry at the equilibrium geometry. As such, it has been studied extensively both from the experimental^90,93–97 and theoretical^34,98–104 points of view. The spectrum is a challenge for computational methods because it is highly resolved, exhibits Herzberg–Teller effects, and contains hot bands due to finite temperature. Although benzene is typically considered to be a rigid molecule, we have recently shown that the anharmonicity affects significantly the intensities of the peaks in the main progression of the spectrum.^25,26 However, our previous work assumed zero temperature, therefore neglecting the weak hot bands present in the experimental spectrum.

Here, we complement our earlier result with the new finite-temperature extended thawed Gaussian method. First, we demonstrate [Fig. 3(a)] the effect of non-zero temperature on the spectrum. Whereas the original, zero-temperature extended thawed Gaussian approximation neglects completely the weak, but non-negligible, hot bands, the finite-temperature approach reproduces all features of the spectrum. The inaccuracy in the frequencies of the peaks is most likely due to the electronic structure method used; we will discuss this later. Nevertheless, Fig. 3(a) clearly shows the difference in the spectra computed without and with finite-temperature effects.

We argue that the benzene absorption spectrum is affected by the anharmonicity of the excited-state potential energy surface. This effect is best demonstrated by the difference in spectra based on two global harmonic models: if the potential energy surface were harmonic, the second-order expansion of the potential energy about any molecular geometry would result in the same spectrum. As shown in Fig. 3(b), in benzene, the adiabatic harmonic method is much more accurate than the vertical; in general, either of the two methods can be more appropriate.^22,23,105 The extended thawed Gaussian approximation outperforms not only the vertical harmonic approach, whose spectrum is completely off, but also the adiabatic harmonic approximation, which fails to produce accurate peak intensities.

Figure 3(c) shows the importance of treating the Herzberg–Teller effect with the extended thawed Gaussian approximation. Since the transition is symmetry-forbidden, i.e., μ(q_eq) = 0, the spectrum computed within the Condon approximation [μ(q) ≈ μ(q_eq)] vanishes, whereas the full, Herzberg–Teller treatment reproduces the experimental spectrum.

In computational chemistry, vibrational scaling factors,¹⁰⁶ which we denote by f, are often used to empirically correct for systematic errors in the vibrational frequencies computed with electronic structure methods. In vibronic spectroscopy, such scaling, applied to ground- and excited-state frequencies, can modify both peak positions and intensities.¹⁰ However, the effect on intensities is often weak; indeed, the adiabatic harmonic spectrum with scaled vibrational frequencies (red, dashed line in Fig. 4) exhibits almost perfect peak positions, but still the same errors in intensities as the adiabatic harmonic spectrum of Fig. 3(b). For comparison—and for comparison only—we show an analogous, “corrected” spectrum computed with the extended thawed Gaussian approximation (blue, solid line in Fig. 4). Since the simple procedure of scaling the vibrational frequencies is not applicable in this case, we scale directly the frequency axis by f, which corrects peak positions but leaves intensities unchanged. The results imply that the subtle anharmonicity effects on spectral intensities, described well with the on-the-fly semiclassical thawed Gaussian method, cannot be captured even with the corrected harmonic potential.

In conclusion, we have presented a new approach to compute vibronic spectra at finite temperature within the framework of the thawed Gaussian approximation. The proposed method describes partially the effect of anharmonicity on the spectrum and at the same time includes all effects treated in the conventional global harmonic approach—mode–mode coupling, non-zero temperature, and Herzberg–Teller contribution to the transition dipole moment. Most importantly, the inclusion of finite temperature comes at no additional computational cost or deterioration in accuracy. Hence, the proposed procedure provides a viable route to systematically improve on global harmonic simulations at any temperature. Finally, this on-the-fly ab initio semiclassical approach to thermo-field dynamics could inspire other quantum or semiclassical “direct dynamics” methods for computing spectra at finite temperatures.

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

The authors would like to thank Maxim Gelin for introducing them to thermo-field dynamics and Lipeng Chen for helpful discussions. 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) is gratefully acknowledged.

Following Suzuki,⁴⁵ let us define

where $|kk̃$ denotes a basis vector of a space obtained as a direct product of “physical” (with basis {|k⟩}) and “fictitious” (with basis${|k̃}$⁠) Hilbert spaces. In general, we use tilde ^̃ to denote an element of (or an operator acting on) the “fictitious” Hilbert space and bar ^¯ (as opposed to bold font used in Ref. 47) for the direct-product space. The physical and fictitious states are related through the conjugation rule⁴⁵ 

which results in

for arbitrary complex numbers u₁ and u₂, states |α⟩ and |α′⟩, and operator Â. Next, the so-called thermal vacuum is defined as

where $ρ^$ is the density operator and acts only on the physical Hilbert space. Then, the correlation function, defined in Eq. (19), can be written as

where $|ϕ¯0=μ^|0¯(β)$ and

The proof goes as follows:

where we used (A3) to go from (A9) to (A10).

To complete the equivalence between Eq. (A5) and Eq. (31) of the main text, we demonstrate that

i.e., the position representation of state $|ϕ¯0$ introduced in this appendix is the function defined in Eq. (23) of the main text. Indeed,

where we again used Eq. (A3) with  being the identity operator.

Here, we derive expressions (32)–(35) for the Gaussian parameters of ρ^1/2(q, q′). First, recall that the matrix element ρ(q, q′) of the thermal density operator in a one-dimensional harmonic oscillator

with mass m and frequency ω is¹⁰⁷ 

Using Eq. (B2), we now derive the expression for the matrix element ρ(q, q′) in a D-dimensional coupled harmonic oscillator,

where m and K are D × D symmetric mass and force-constant matrices, respectively, and q_eq is the equilibrium position at which the potential energy has its minimum.

The first step consists in transforming q to the mass-scaled normal-mode coordinates Q = O^T · m^1/2 · (q − q_eq), where O is the orthogonal matrix diagonalizing the mass-scaled force-constant matrix, i.e., $OT⋅m−1/2⋅K⋅m−1/2⋅O=Ωdiag2$ with a real diagonal matrix Ω_diag. This leads to the uncoupled Hamiltonian,

where $Ĥi=12P^i2+12ωi2Qi^2$ and ω_i are the diagonal elements of Ω_diag. The density matrix element of this uncoupled D-dimensional harmonic oscillator is

where

are 2D × 2D matrices composed of D × D diagonal sub-matrices.

Transformation back to the original coordinates q yields

with

To find ρ^1/2(q, q′), we rewrite it as

and hence, the initial value Ā₀ needed for the extended thawed Gaussian propagation is given by Ā_β/2, where Ā_β is defined by Eqs. (B16), (B18), and (B19), which proves Eqs. (32)–(34). To find the initial value $γ¯0$⁠, we can avoid computing the traces in Eq. (B20) explicitly and instead recognize that ρ^1/2(q, q′) must be normalized as

Therefore, $γ¯0$ can be computed from Ā₀ in analogy to Eq. (14) for the wavepacket (3), i.e.,

where we used that the initial width matrix Ā₀ is Ā_β/2 and that

(in particular, the determinant is independent of ℏ and temperature) because

In this derivation, we used the relation¹⁰⁸ 

valid for arbitrary invertible matrices A and B, the relation