Sampling strategies for the Herman–Kluk propagator of the wavefunction

When the semiclassical Herman–Kluk propagator is used for evaluating quantum-mechanical observables or time-correlation functions, the initial conditions for the guiding trajectories are typically sampled from the Husimi density. Here, we employ this propagator to evolve the wavefunction itself. We investigate two grid-free strategies for the initial sampling of the Herman–Kluk propagator applied to the wavefunction and validate the resulting time-dependent wavefunctions evolved in harmonic and anharmonic potentials. In particular, we consider Monte Carlo quadratures based either on the initial Husimi density or on its square root as possible and most natural sampling densities. We prove analytical convergence error estimates and validate them with numerical experiments on the harmonic oscillator and on a series of Morse potentials with increasing anharmonicity. In all cases, sampling from the square root of Husimi density leads to faster convergence of the wavefunction.

1 Introduction

Semiclassical initial value representation techniques [1, 2] have evolved into useful tools for calculations of the dynamics of atoms and molecules [3]. Frozen Gaussians and their superposition were introduced by Heller in 1981 [4] as an extension to the thawed Gaussian approximation [5] in order to capture non-linear spreading of the wavepackets. Herman and Kluk justified the frozen-Gaussian ansatz and introduced an improved approximation [6–8], now known as the Herman–Kluk propagator, which contains an additional prefactor that rigorously compensates for the fixed width of these frozen Gaussians and ensures unitarity of the time evolution in the stationary-phase limit [7]. The Herman–Kluk approximation keeps the trajectories of the individual Gaussians uncoupled, which sets it apart from the more accurate, but computationally more demanding approaches like the coupled coherent states [9] or the variational multi-configurational Gaussian method [10].

In the spirit of the semiclassical initial value representation, the Herman–Kluk propagator avoids the root search problem [2, 11]. Because of its high accuracy, this propagator belongs among the most successful semiclassical approximations [12–18] and has been derived in many different ways [9, 19–22]. There are observables, such as low-resolution vibronic spectra in mildly anharmonic systems, which can be well described by the more efficient thawed Gaussian approximation [23, 24] and other single-trajectory methods [25–30]. In chaotic and other systems, however, increased anharmonicity leads to wavepacket splitting and non-trivial interference effects. In such situations, the single-trajectory methods break down, whereas the Herman–Kluk propagator often remains accurate and, at the least, provides a qualitatively correct insight (see, e.g., Figure 1). However, in its basic form, the Herman–Kluk propagator cannot describe non-adiabatic dynamics and quantum tunneling. To overcome this limitation, various extensions have been developed: For example, surface hopping Herman–Kluk initial value representation [31, 32] captures non-adiabatic dynamics, whereas higher-order semiclassical corrections to the Herman–Kluk propagator [33, 34] incorporate nuclear tunneling.

FIGURE 1

FIGURE 1. Upper panel: Spectra of a Morse potential evaluated using the exact quantum dynamics, Herman–Kluk (HK) propagator, and thawed Gaussian approximation (TGA). Both approximations yield accurate results. Lower panel: Position density at time t ≈ 392 fs propagated in the same Morse potential. In contrast to the Herman–Kluk propagator, the thawed Gaussian approximation does not capture interference between faster and slower components of the wavepacket. For more details, see the last paragraph of Section 4.

The multi-trajectory nature of the Herman–Kluk propagator is not only its advantage, but also its bottleneck. In particular, converged computations might require an extremely large number of trajectories. For this reason, several groups designed various methods, whose goal is reducing the number of trajectories required by the Herman–Kluk propagator. These include, e.g., Filinov filtering [13, 35–37], time-averaging [38–41], semiclassical interaction picture [42, 43], multiple-coherent states [44], hybrid dynamics [45–47], mixed quantum-semiclassical dynamics [48–50], and many others, which achieve the reduction in the number of trajectories by applying one of several possible further approximations.

Yet, the number of guiding trajectories can already be significantly reduced simply by choosing the sampling density of the initial conditions wisely. Although the acceleration of the convergence may be smaller than with the previously mentioned approximate methods, the advantage of this “purely numerical” approach based on improved exact sampling of the unaltered Herman–Kluk initial value representation is that the converged results agree exactly with the converged results of the original Herman–Kluk propagator. For this reason, an early numerical study [8] of the Herman–Kluk wavefunction employed the square root of the Husimi distribution as the sampling density of the initial state, whereas most calculations of observables, time-correlation functions, and wavepacket autocorrelation functions, i.e., quantities quadratic in the wavefunction, correctly employ sampling from the Husimi density [2, 51]. In contrast to observables and correlation functions, the Herman–Kluk wavefunction itself has not been studied much since the early papers [6, 8] and rigorous numerical analysis has been presented only recently [52, 53]. However, the choice of the optimal sampling density for the Herman–Kluk wavefunction itself has not been analyzed in detail.

The goal of this work is, therefore, to analyze the convergence of the Herman–Kluk wavefunction for different initial sampling strategies and to understand the convergence error as a function of time. Specifically, we investigate two mesh-free discretization approaches for the initial sampling—first analytically and then numerically on the examples of harmonic and Morse oscillators with increasing anharmonicity. In a follow-up paper, an analogous detailed analysis of the convergence of the norm, energy, and other expectation values will be presented.

The remainder of this paper is organized as follows. In the next section we briefly introduce the Herman–Kluk propagator and its components necessary for numerical computations. We define the initial sampling densities and specify the algorithm used for numerical experiments. In the main Section 3, we analyze the errors at the initial and final times due to the phase-space discretization. In particular, we prove that in the harmonic oscillator this error is a periodic function of time. Our numerical experiments in Section 4 confirm the theoretical error estimates and provide further insights into the anharmonic evolution generated by Morse potentials, which are not accessible to explicit analytical calculations.

2 Discretising the Herman–Kluk propagator

2.1 Herman–Kluk propagator

Evolution of a quantum state |ψ(t)⟩ is governed by the time-dependent Schrödinger equation

$i\hslash \frac{d}{dt}|\psi \left(t\right)⟩=\widehat{H}|\psi \left(t\right)⟩,|\psi \left(0\right)⟩=|{\psi }_0⟩(1)$

where ℏ is the reduced Planck constant and $\widehat{H}$ is the Hamiltonian operator. Here, we assume the Hamiltonian to be the time-independent operator

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

where m is the mass (after mass-scaling coordinates), and q and p are D-dimensional position and momentum vectors. Under the usual regularity and growth assumptions on the potential energy function V, the spectral-theorem [54] provides for all times $t\in \mathbb{R}$ a well-defined unitary propagator

${\widehat{U}}_t:=\exp \left(-i\widehat{H}t/\hslash \right),(3)$

in terms of which the solution of (1) can be expressed as

$|\psi \left(t\right)⟩={\widehat{U}}_t|{\psi }_0\rangle (4)$

for all square integrable initial data ${\psi }_0(x)=\langle x|{\psi }_0\rangle \in L^2({\mathbb{R}}^D)$ .

Solving the Schrödinger equation numerically is notoriously difficult for various reasons. With respect to the atomic scale, the nuclear mass is rather large, so that the presence of the factor m ⁻¹ in the kinetic energy operator induces highly oscillatory motion both with respect to time and space. More importantly, for most molecular systems the dimension D of the configuration space is so large that grid-based integration methods are very expensive if not infeasible. Therefore, one often resorts to mesh-free discretization methods or semiclassical approximations, both of which alleviate this “curse of dimensionality” at least partially.

The semiclassical Herman–Kluk propagator utilizes frozen Gaussian functions

$\begin{array}{c}{g}_z^{\gamma }\left(x\right)={\left(\frac{\det \gamma }{{\pi }^D{\hslash }^D}\right)}^{1/4}\hfill \\ \times \exp \left\{-\left[{\left(x-q\right)}^T\cdot \gamma \cdot \left(x-q\right)/2+i{p}^T\cdot \left(x-q\right)\right]/\hslash \right\}\hfill \end{array}(5)$

with a centre at the phase-space point $z=(q,p)\in {\mathbb{R}}^{2D}$ and with a fixed, real, symmetric, positive-definite width-matrix $\gamma \in {\mathbb{R}}^{D\times D}$ . The frozen Gaussians $\{|g_z^{\gamma }\rangle :z\in {\mathbb{R}}^{2D}\}$ form an over-complete subset of the Hilbert space of square integrable functions. They have the striking property that any $|\psi \rangle \in L^2({\mathbb{R}}^D)$ can be decomposed as [55]

$|\psi ⟩={\int }_{{\mathbb{R}}^{2D}}\langle g_z^{\gamma }|\psi \rangle |g_z^{\gamma }\rangle d\nu ,(6)$

with the scaled phase-space measure dν = dz/(2πℏ) ^D . Using this decomposition for our solution of the time-dependent Schrödinger equation, we obtain

${\widehat{U}}_t|{\psi }_0⟩={\int }_{{\mathbb{R}}^{2D}}\langle g_z^{\gamma }|{\psi }_0\rangle {\widehat{U}}_t|g_z^{\gamma }\rangle d\nu .(7)$

Approximating the exact propagator ${\widehat{U}}_t$ with the Herman–Kluk propagator

${\widehat{U}}_t^{\text{HK}}:={\int }_{{\mathbb{R}}^{2D}}R\left(t,z\right)e^{iS\left(t,z\right)/\hslash }|g_{z\left(t\right)}^{\gamma }\rangle \langle g_z^{\gamma }|d\nu (8)$

yields the Herman–Kluk wave function

$\begin{array}{ll}\hfill |{\psi }^{\text{HK}}\left(t\right)\rangle & :={\widehat{U}}_t^{\text{HK}}|{\psi }_0\rangle \hfill \\ \hfill & ={\int }_{{\mathbb{R}}^{2D}}R\left(t,z\right)e^{iS\left(t,z\right)/\hslash }\langle g_z^{\gamma }|{\psi }_0\rangle |g_{z\left(t\right)}^{\gamma }\rangle d\nu .\hfill \end{array}(9)$

Here, z(t) = (q(t), p(t)) is the solution to the underlying classical Hamiltonian system

$\dot{w}=J\cdot \nabla h\left(w\right),w\left(0\right)=z,(10)$

for the Hamiltonian function h(q, p) = T(p) + V(q), where

$J=\left(\begin{array}{cc}\hfill 0\hfill & \hfill {\mathrm{I}\mathrm{d}}_D\hfill \\ \hfill -{\mathrm{I}\mathrm{d}}_D\hfill & \hfill 0\hfill \end{array}\right)\in {\mathbb{R}}^{2D\times 2D}(11)$

is the standard symplectic matrix. S denotes the classical action integral

$S\left(t,z\right)={\int }_0^t\left[\frac{d}{d\tau }q\left(\tau \right)\cdot p\left(\tau \right)-h\left(z\left(\tau \right)\right)\right]d\tau (12)$

which solves the initial value problem

$\dot{S}\left(t,z\right)=T\left(p\left(t\right)\right)-V\left(q\left(t\right)\right),S\left(0,z\right)=0,(13)$

for all $z=(q,p)\in {\mathbb{R}}^{2D}$ . In (9), the Herman–Kluk prefactor

$R\left(t,z\right)=2^{-D/2}\cdot \det {\left(M_{qq}+{\gamma }^{-1}\cdot M_{pp}\cdot \gamma -i{M}_{qp}\cdot \gamma +i{\gamma }^{-1}\cdot M_{pq}\right)}^{1/2}(14)$

depends on the matrices

$M_{\alpha \beta }=\partial {\alpha }_t/\partial {\beta }_0\in {\mathbb{R}}^{D\times D}\left(\alpha ,\beta \in \left\{q,p\right\}\right),(15)$

i.e., the four D × D block components of the 2D × 2D stability matrix M. Stability matrix, defined as the Jacobian M(t) = ∂z(t)/∂z of the flow map, is the solution to the variational equation

$\dot{M}\left(t\right)=J\cdot \mathrm{H}\mathrm{e}\mathrm{s}\mathrm{s}h\left(z\left(t\right)\right)\cdot M\left(t\right),M\left(0\right)={\mathrm{I}\mathrm{d}}_{2D},(16)$

with Hess h the Hessian of the Hamiltonian function h.

The Herman–Kluk approximation (9) has been mathematically justified in different works [53, 56–58]. It has been shown that the exactly evolved quantum state is approximated by the Herman–Kluk state (9) with an error of the order of ℏ. More precisely,

$\underset{t\in \left[0,T\right]}{\sup }{|||\psi \left(t\right)\rangle -|{\psi }^{\text{HK}}\left(t\right)\rangle ||}_{L^2}\le C\left(T\right)\hslash ,(17)$

for all initial data |ψ ₀〉 of norm one, where T > 0 is a fixed time and C(T) > 0 is a constant independent of ℏ and independent of |ψ ₀〉. If the potential is at most quadratic, then the approximation is exact [53]. For the expectation value of an observable $\widehat{A}$ , the error of the Herman–Kluk approximation can be pessimistically estimated as

$\langle \psi \left(t\right)|\widehat{A}|\psi \left(t\right)\rangle =\langle {\psi }^{\text{HK}}\left(t\right)|\widehat{A}|{\psi }^{\text{HK}}\left(t\right)\rangle +\mathcal{O}\left(\hslash \right),(18)$

by using the triangle inequality. In particular, this gives an upper bound for the error of the squared norm (if $\widehat{A}=\widehat{I}$ ) or energy (if $\widehat{A}=\widehat{H}$ ). However, this coarse estimate is potentially not sharp, since it cannot account for error cancellation due to oscillation.

2.2 Discretisation

Evolving the wavefunction (9) with the Herman–Kluk propagator requires evaluating an integral over the phase space ${\mathbb{R}}^{2D}$ and the overlap of the initial state with a frozen Gaussian. Furthermore, the algorithm needs to propagate the trajectories z(t) = (q(t), p(t)), the classical action S(t, z), and the Herman–Kluk prefactor R(t, z) according to Hamilton’s equations of motion for all chosen quadrature points $z\in {\mathbb{R}}^{2D}$ . The latter can be achieved using symplectic integration methods to preserve also the symplectic structure of the classical Hamiltonian system. Due to the curse of dimensionality, for high D the integral on ${\mathbb{R}}^{2D}$ must be evaluated using mesh-free discretization, such as Monte Carlo methods [59].

To evaluate the Herman–Kluk wavefunction (9) by Monte Carlo sampling, we rewrite the integral from Eq. 9 as

$|{\psi }^{\text{HK}}\left(t\right)\rangle ={⟨r\left(z\right)|\varphi \left(t,z\right)⟩⟩}_{\rho \left(z\right)},(19)$

where we introduced the notation

$\begin{array}{cc}\hfill {⟨|\psi \left(t,z\right)⟩⟩}_{\rho \left(z\right)}& :={\int }_{{\mathbb{R}}^{2D}}|\psi \left(t,z\right)⟩d\mu \hfill \\ \hfill & ={\int }_{{\mathbb{R}}^{2D}}|\psi \left(t,z\right)\rangle \rho \left(z\right)d\nu \hfill \end{array}(20)$

for the phase-space average of |ψ(t, z)⟩ with respect to a probability measure μ with density ρ(z) = dμ/dν and defined a time-dependent state

$|\varphi \left(t,z\right)⟩:=R\left(t,z\right)e^{iS\left(t,z\right)/\hslash }|g_{z\left(t\right)}^{\gamma }\rangle ,(21)$

i.e., a propagated frozen Gaussian multiplied by the Herman–Kluk and phase prefactors, and a time-independent function

$r\left(z\right):=\langle g_z^{\gamma }|{\psi }_0\rangle \rho {\left(z\right)}^{-1}.(22)$

The Monte Carlo estimator is then given by

$|{\psi }_N\left(t\right)\rangle =\frac{1}{N}\sum _{j=1}^N|\psi (t,z_j)\rangle ,(23)$

where ψ(t, z _j ) is the state obtained from initial condition z _j and z ₁, z ₂, … , z _N are sampled from probability density ρ(z).

As for any importance sampling, there are infinitely many ways to decompose the time-independent part of the phase-space integrand in Eq. 9 into the product $\langle g_z^{\gamma }|{\psi }_0\rangle =r(z)\rho (z)$ of a prefactor r with a normalized sampling density ρ. If one computes observables and correlation functions, which are quadratic in the initial state, ρ(z) is typically taken to be the Husimi probability density

${\rho }_{\mathrm{H}}\left(z\right):=|\langle g_z^{\gamma }|{\psi }_0\rangle {|}^2(24)$

of the initial state |ψ ₀〉. For Husimi distribution, the prefactor r(z) becomes

$r_{\mathrm{H}}\left(z\right)={\langle {\psi }_0|g_z^{\gamma }\rangle }^{-1}.(25)$

However, since the wavefunction evolved with the Herman–Kluk propagator is first- and not second-order in |ψ ₀〉, it is natural to also test the square root of the Husimi density and to consider the probability density

$\tilde{\rho }\left(z\right):=\frac{|\langle g_z^{\gamma }|{\psi }_0\rangle |}{{\int }_{{\mathbb{R}}^{2D}}|\langle g_w^{\gamma }|{\psi }_0\rangle |d\nu }(26)$

with a bounded prefactor

$\tilde{r}\left(z\right):=\langle g_z^{\gamma }|{\psi }_0\rangle \tilde{\rho }{\left(z\right)}^{-1}.(27)$

Indeed, in their early paper [8], Kluk et al. used this “square-root” approach for computing the Herman–Kluk wavefunction and mentioned that it was “especially appealing because it [was] a well defined, non-arbitrary way of choosing the basis […].” In the following, we shall show that the square-root approach is indeed optimal, but the number of samples still suffers from an exponential dependence on the dimension.

The two probability densities ρ _H(z) and $\tilde{\rho }(z)$ can now be used to compute the phase-space integral by Monte Carlo integration. To sum up, we have two cases, in which we evaluate |ψ ^HK(t)〉 either as

$|{\psi }^{\text{HK}}\left(t\right)\rangle ={\langle r_{\mathrm{H}}\left(z\right)|\varphi \left(t,z\right)\rangle \rangle }_{{\rho }_{\mathrm{H}}\left(z\right)}\mathrm{(Case\; H)}$

or as

$|{\psi }^{\text{HK}}\left(t\right)\rangle ={\langle \tilde{r}\left(z\right)|\varphi \left(t,z\right)\rangle \rangle }_{\tilde{\rho }\left(z\right)}.\mathrm{(Case\; sqrt\; ‒\; H)}$

In general, the integral over ${\mathbb{R}}^D$ defining the overlap of the initial wavefunction with a Gaussian has to be computed by numerical quadrature. However, for important specific cases analytical formulas are available. If the initial wavefunction is a Gaussian wavepacket $|{\psi }_0\rangle =|g_{z_0}^{\gamma }\rangle$ centred at some phase-space point $z_0=(q_0,p_0)\in {\mathbb{R}}^{2D}$ , then

$\begin{array}{cc}\hfill {\langle g_z^{\gamma }|{\psi }_0\rangle }_{L^2\left({\mathbb{R}}^D\right)}& =\exp \left[-\frac{1}{4\hslash }{\left(z-z_0\right)}^T\cdot {\mathrm{\Sigma }}_0\cdot \left(z-z_0\right)\right]\hfill \\ \hfill & \cdot \exp \left[\frac{i}{2\hslash }{\left(p+p_0\right)}^T\cdot \left(q-q_0\right)\right],\hfill \end{array}(28)$

where ${\mathrm{\Sigma }}_0:\; =\left(\begin{array}{cc}\hfill \gamma \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\gamma }^{-1}\hfill \end{array}\right)$ is the matrix containing the width parameters of the initial coherent state. Then, the Husimi density is given by

${\rho }_{\mathrm{H}}\left(z\right)=\exp \left[-\frac{1}{2\hslash }{\left(z-z_0\right)}^T\cdot {\mathrm{\Sigma }}_0\cdot \left(z-z_0\right)\right],(29)$

whereas the second approach provides the density

$\tilde{\rho }\left(z\right)=2^{-D}\exp \left[-\frac{1}{4\hslash }{\left(z-z_0\right)}^T\cdot {\mathrm{\Sigma }}_0\cdot \left(z-z_0\right)\right](30)$

and prefactor

$\tilde{r}\left(z\right)=2^D\exp \left[\frac{i}{2\hslash }{\left(p+p_0\right)}^T\cdot \left(q-q_0\right)\right].(31)$

2.3 Summary of the numerical algorithm

Taking into account all the previous considerations, we slightly extend the natural numerical algorithm (described, e.g., Section 4 of Ref. [52]) for finding the Herman–Kluk approximation to the wavefunction at time t.

Algorithm 1:. (Herman–Kluk propagation)

1. Draw independent samples $z_1,\dots ,z_N\in {\mathbb{R}}^{2D}$ from a distribution with density ρ _H(z) or $\tilde{\rho }(z)$ given by Eqs. 29 and 30.

2. For all j ∈ {0, 1, … , N}:

2.1 Set initial values z(0) = z _j , M(0) = Id_2D and S(0) = 0.

2.2 Compute approximate solutions to Eqs. 10, 13 and 16 up to time t with a symplectic integration method [60] based on the Störmer–Verlet scheme [61, 62].

2.3 Compute the Herman–Kluk prefactor R(t, z _j ) from M(t) while choosing the correct branch of the complex square root [63], which guarantees continuity of R(t, z _j ) as a function of t.

3. Calculate |ψ _N (t)〉 by means of formula (23) with |ψ(t, z _j )〉 replaced with either

$|{\psi }_{\mathrm{H}}(t,z_j)⟩=r_{\mathrm{H}}(z_j)R(t,z_j)e^{iS\left(t,z_j\right)/\hslash }|g_{z_j\left(t\right)}^{\gamma }⟩(32)$

$\begin{array}{ll}\hfill \text{or}& \hfill \\ \hfill & |\tilde{\psi }(t,z_j)\rangle =\tilde{r}(z_j)R(t,z_j)e^{iS\left(t,z_j\right)/\hslash }|g_{z_j\left(t\right)}^{\gamma }\rangle ,\hfill \end{array}(33)$

where r _H(z _j ) and $\tilde{r}(z_j)$ are given by Eqs. 25, 28 and 31.

We note that Eqs. 10, 13 and 16 can be evaluated simultaneously with a single numerical integrator. To increase the accuracy of the time integration one can use higher-order composition methods [64, 65]. They increase the order of the time integrator, but also its numerical cost. A higher-order time integrator is not necessarily useful since the phase-space error occurring from the Monte Carlo quadrature usually dominates the time integration error.

3 Phase space error analysis

Algorithm 1 relies on the discretisation of the phase-space integrals and of the system of ordinary differential equations. Here we focus only on the phase-space discretisation error. A similar, but less formal phase-space error analysis was applied to various algorithms for computing fidelity and classical correlation functions [66–69].

3.1 Moments of the integrand

To assess the accuracy of the Monte Carlo estimator (23), we examine the moments of the two different integrands: |ψ _H(t)〉 and $|\tilde{\psi }(t)\rangle$ . First, we observe that for any s > 0,

$\mathbb{E}\left[\Vert \tilde{\psi }\left(t\right){\Vert }^s\right]={\int }_{{\mathbb{R}}^{2D}}|R\left(t,z\right){|}^s|\tilde{r}\left(z\right){|}^{s}d\tilde{\mu }\left(z\right)<\infty ,(34)$

since $\tilde{r}(z)$ and the Herman–Kluk prefactor R(t, z) are both bounded functions. For a discussion of the boundedness of R(t, z), see Section 1 of the Supplementary Material. In contrast,

$\begin{array}{cc}\hfill \mathbb{E}\left[\Vert {\psi }_{\mathrm{H}}\left(t\right){\Vert }^s\right]& ={\int }_{{\mathbb{R}}^{2D}}|R\left(t,z\right){|}^s|r_{\mathrm{H}}\left(z\right){|}^{s}d{\mu }_{\mathrm{H}}\left(z\right)\hfill \\ \hfill & ={\int }_{{\mathbb{R}}^{2D}}|R\left(t,z\right){|}^s|\langle g_z^{\gamma }|{\psi }_0\rangle {|}^{2-s}d\nu \hfill \end{array}(35)$

ceases to be finite for s ≥ 2. However, both integrands have a finite first moment, so that the strong law of large numbers (Theorem 2.4.1 in Ref. [70]) provides convergence of the estimator,

$|{\psi }_N\left(t\right)\rangle \to |{\psi }^{\text{HK}}\left(t\right)\rangle \text{as}N\to \infty ,(36)$

with probability one. Divergence of the second moment for (Case H) does not violate convergence of the estimator, but results in a slightly worse convergence rate than the one for (Case sqrt-H). Numerical results in Section 4 confirm this expectation. This shows that the initial sampling density has to be chosen carefully.

3.2 Mean squared error

For (Case sqrt-H) the second moment is finite, so that the mean squared error of the Monte Carlo estimator (23) is well-defined and satisfies

$\mathbb{E}\left[\Vert |{\psi }_N\left(t\right)\rangle -|{\psi }^{\text{HK}}\left(t\right)\rangle {\Vert }^2\right]=\frac{\mathbb{V}\left[|\tilde{\psi }\left(t\right)\rangle \right]}{N},(37)$

where the expectation value and the variance are with respect to the density $\tilde{\rho }(z)$ . Moreover (see also [52]),

$\begin{array}{ll}\hfill \mathbb{V}\left[|\tilde{\psi }\left(t\right)\rangle \right]& ={\int }_{{\mathbb{R}}^{2D}}\Vert |\tilde{\psi }\left(t\right)\rangle -\mathbb{E}\left[|\tilde{\psi }\left(t\right)\rangle \right]{\Vert }^{2}d\tilde{\mu }\left(z\right)\hfill \\ \hfill & ={\int }_{{\mathbb{R}}^{2D}}|R\left(t,z\right){|}^2|\tilde{r}\left(z\right){|}^{2}d\tilde{\mu }\left(z\right)-\Vert |{\psi }^{\text{HK}}\left(t\right)\rangle {\Vert }^2.\hfill \end{array}(38)$

In the special case of an initial Gaussian initial wavepacket $|{\psi }_0\rangle =|g_{z_0}^{\gamma }\rangle$ , this simplifies to

$\mathbb{V}\left[|\tilde{\psi }\left(t\right)\rangle \right]=4^D{\int }_{{\mathbb{R}}^{2D}}|R\left(t,z\right){|}^{2}d\tilde{\mu }\left(z\right)-\Vert |{\psi }^{\text{HK}}\left(t\right)\rangle {\Vert }^2,(39)$

and at initial time t = 0 we obtain

$\mathbb{V}\left[|\tilde{\psi }\left(0\right)\rangle \right]=4^D-1,(40)$

since the Herman–Kluk prefactor satisfies R(0, z) = 1.

For a more general assessment of the variance, numerical experiments in Section 4 consider the error between the approximations with N and 2N samples. By the linearity of the expectation value and the triangle inequality, this error can be estimated by

$\begin{array}{c}\mathbb{E}\left[\Vert |{\psi }_N\left(t\right)\rangle -|{\psi }_{2N}\left(t\right)\rangle {\Vert }^2\right]\hfill \\ \le \mathbb{E}\left[\Vert |{\psi }_N\left(t\right)\rangle -|{\psi }^{\text{HK}}\left(t\right)\rangle {\Vert }^2\right]+\mathbb{E}\left[\Vert |{\psi }^{\text{HK}}\left(t\right)\rangle -|{\psi }_{2N}\left(t\right)\rangle {\Vert }^2\right]\hfill \\ =\frac{3}{2}\frac{\mathbb{V}\left[|\psi \left(t\right)\rangle \right]}{N}.\hfill \end{array}(41)$

Note that Eq. 37 gives only the expected convergence error, i.e., the convergence error averaged over infinitely many independent simulations, each using N trajectories. The actual convergence error for any specific simulation with N trajectories may deviate from this analytical estimate substantially due to statistical noise. Nevertheless, in Section 2 of the Supplementary Material, we explain how Eq. 37 can also provide a rigorous lower bound and asymptotic estimate of the number of trajectories needed for convergence of a single simulation.

3.3 Other sampling densities

For the special case of an initial Gaussian wavepacket $|{\psi }_0\rangle =|g_{z_0}^{\gamma }\rangle$ , the two proposed sampling densities ρ _H(z) and $\tilde{\rho }(z)$ belong to a family of normal distributions with probability density functions

${\rho }_a\left(z\right)={\left(\frac{2}{a}\right)}^D\exp \left[-\frac{1}{a\hslash }{\left(z-z_0\right)}^T\cdot {\mathrm{\Sigma }}_0\cdot \left(z-z_0\right)\right](42)$

with a ≥ 2. In the spirit of importance sampling, the Herman–Kluk wavefunction can accordingly be written as a phase-space average

$|{\psi }^{\text{HK}}\left(t\right)\rangle ={\langle |{\psi }_a\left(t,z\right)\rangle \rangle }_{{\rho }_a\left(z\right)}.(43)$

At time t = 0, the norm of the integrand satisfies

$\Vert |{\psi }_a\left(0,z\right)\rangle \Vert ={\left(\frac{a}{2}\right)}^D\exp \left[-\left(\frac{a-4}{4a\hslash }\right){\left(z-z_0\right)}^T\cdot {\mathrm{\Sigma }}_0\cdot \left(z-z_0\right)\right],(44)$

which implies for the variance

$\begin{array}{ll}\hfill \mathbb{V}\left[|{\psi }_a\left(0\right)\rangle \right]& ={\left(\frac{a}{4\pi \hslash }\right)}^D{\int }_{{\mathbb{R}}^{2d}}\exp \left[-\left(\frac{a-2}{2a\hslash }\right)|z{|}^2\right]dz-1\hfill \\ \hfill & ={\left(\frac{a^2}{2\left(a-2\right)}\right)}^D-1.\hfill \end{array}(45)$

For a > 2, the function a↦a ²/(2a − 4) attains its minimum at a = 4, which corresponds to the sampling density $\tilde{\rho }$ . In other words, (Case sqrt-H) is optimal as far as the mean squared error is concerned. Nevertheless, even this optimal sampling results in an unfavorable exponential growth with the number of dimensions D.

3.4 Harmonic motion

The one-dimensional harmonic oscillator

$\widehat{H}=-\frac{{\hslash }^2}{2m}\frac{d^2}{d{x}^2}+\frac{m{\omega }^2}{2}{x}^2(46)$

is one of the rare examples for which explicit expressions for the solution of the Schrödinger equation and for the Herman–Kluk prefactor exist. For an initial Gaussian wavepacket $|{\psi }_0\rangle =|g_{z_0}^{\gamma }\rangle$ with position and momentum $q_0,p_0\in \mathbb{R}$ , the exact wavefunction is given by [71]

${\psi }_{\text{ex}}\left(x,t\right)=\exp \left\{\frac{i}{\hslash }\left[\frac{{\alpha }_t}{2}{\left(x-q_t\right)}^2+p_t\left(x-q_t\right)+{\beta }_t\right]\right\},(47)$

where

${\alpha }_t=b\frac{{\alpha }_0\cos \left(\omega t\right)-b\sin \left(\omega t\right)}{b\cos \left(\omega t\right)+{\alpha }_0\sin \left(\omega t\right)},(48)$

$q_t=q_0\cos \left(\omega t\right)+\frac{p_0}{b}\sin \left(\omega t\right),(49)$

$p_t=p_0\cos \left(\omega t\right)-b{q}_0\sin \left(\omega t\right)\text{\hspace{0.17em}and\hspace{0.17em}}(50)$

${\beta }_t={\beta }_0+\frac{1}{2}\left[q_t{p}_t-q_0{p}_0+i\hslash \ln \left(\frac{z_t}{b}\right)\right](51)$

with the abbreviations z _t = b cos(ωt) + α ₀ sin(ωt), α ₀ = iγ, and b = mω. Here, β ₀ includes the normalization constant for the wavefunction at time t = 0. The classical action is given by

$S\left(t,q_0,p_0\right)=\frac{1}{2}\sin \left(\omega t\right)\left[\left(\frac{p_0^2}{b}-b{q}_0^2\right)\cos \left(\omega t\right)-q_0{p}_0\sin \left(\omega t\right)\right].(52)$

The four components of stability matrix M can be obtained, from their definition (15), by differentiating expressions for q _t and p _t with respect to q ₀ and p ₀, namely,

$M\left(t\right)=\left(\begin{array}{cc}\hfill \cos \left(\omega t\right)\hfill & \hfill b^{-1}\sin \left(\omega t\right)\hfill \\ \hfill -b\sin \left(\omega t\right)\hfill & \hfill \cos \left(\omega t\right)\hfill \end{array}\right).(53)$

The Herman–Kluk prefactor satisfies

$R\left(t\right)={\left\{\frac{1}{2}\left[2\cos \omega t-i\sin \left(\omega t\right)\left(\frac{\gamma }{b}+\frac{b}{\gamma }\right)\right]\right\}}^{1/2},(54)$

so that the variance (38) can be written as

$\mathbb{V}\left[|\tilde{\psi }\left(t\right)\rangle \right]=2{\left[4{\cos }^2\left(\omega t\right)+{\left(\frac{\gamma }{b}+\frac{b}{\gamma }\right)}^2{\sin }^2\left(\omega t\right)\right]}^{1/2}-1.(55)$

This implies that for (Case sqrt-H) applied to a harmonic oscillator, the mean squared error of our Monte Carlo estimator oscillates with a period of π/ω between

$3\le \mathbb{V}\left[|\tilde{\psi }\left(t\right)\rangle \right]\le 2\left(\frac{\gamma }{m\omega }+\frac{m\omega }{\gamma }\right)-1.(56)$

4 Numerical examples

In this section, we complement our previous theoretical results with numerical examples. We start by examining how the performance of Algorithm 1 depends on the method to discretise the phase space. We explore the time dependence of the variance of the Monte Carlo estimator in one dimension in a harmonic oscillator as well as in a series of increasingly anharmonic Morse potentials. The Monte Carlo integration is tested by averaging over N independent, identically distributed samples of initial conditions. We approximate its convergence rate by assuming a power law

$F\left(N\right)=c{N}^{-s}(57)$

dependence of the mean statistical error on the number of samples N. The prefactor c and order s of convergence are determined by the linear fit (in the least-squares sense) of the logarithm of Eq. 57 to the dependence of the logarithm of the statistical error on the logarithm of N.

Throughout our numerical examples, we work in atomic units (ℏ = 1), mass-scaled coordinates and with an initial state that is a Gaussian wavepacket with phase-space centre $z_0\in {\mathbb{R}}^{2D}$ .

4.1 Initial phase-space error

We start by considering a spherical initial Gaussian wavepacket $|{\psi }_{\text{ex}}(0)\rangle =|g_{z_0}^{\gamma }\rangle$ with a width parameter γ = 2Id _D in one and four dimensions. For D = 1, it is centred at z ₀ = (−1, 0) and for D = 4, at z ₀ = −(1, 1, 1, 1, 0, 0, 0, 0). Figure 2 shows the L ²-distance

$\Vert |{\psi }_N\left(0\right)\rangle -|{\psi }_{\text{ex}}\left(0\right)\rangle {\Vert }_{L^2}(58)$

of the Monte Carlo estimators for (Case H) and (Case sqrt-H) from the exact wavefunction at initial time as a function of the number of Monte Carlo quadrature points. We can immediately see that the analytical prediction (40) of the mean squared error (37) of (Case sqrt-H) is fulfilled. For (Case H), the curve-fitting approximation (57) provided (c, s) = (2.56, −0.48) for D = 1 and (c, s) = (19.3, −0.36) for D = 4. It shows that both cases converge to the correct result and that (Case sqrt-H) performs slightly better than (Case H). Additionally, Figure 3 displays the L ²-distance of the estimators for both cases from the exact wavefunction at initial time as a function of the dimension D. Each wavefunction was approximated with N ≈ 8 ⋅ 10⁵ trajectories. The analytical prediction (40) for (Case sqrt-H) is realized. Moreover, the error for (Case H) increases faster with D.

FIGURE 2

FIGURE 2. Sampling error of the initial wavefunction in one (upper panel) and four (lower panel) dimensions as a function of the number N of Monte Carlo points. The plot displays the error for sampling from the Husimi density (Case H) and its approximated convergence (marked lines) as well as the error for sampling from the square root of the Husimi density (Case sqrt-H) and its analytical error estimation (dotted line).

FIGURE 3

FIGURE 3. Dependence of the sampling error of the initial wavefunction on dimension D for N = 100 ⋅ 2¹³ ≈ 8 ⋅ 10⁵ points sampled from either the Husimi density (Case H) or its square root (Case sqrt-H). The analytical error estimate for the latter sampling is shown by the dotted line.

4.2 Harmonic potential

To analyze the effect of dynamics on the convergence, let us consider a harmonic potential V(x) = x ²/2 in one dimension and explore the dynamics for one full oscillation period. The initial Gaussian wavepacket is still localized in z ₀ = (−1, 0) with width γ = 2. In Figure 4, we compare the exact wavefunction (47) with the numerical realization of the Herman–Kluk propagator which uses the exact solution to q(t), p(t), S(t, z) and R(t, z) stated in Eqs 49, 50, 52, 54. Since the Herman–Kluk propagator is exact for harmonic motion and since we supply exact classical trajectories, the observed numerical error is only due to the Monte Carlo integration. The upper panel of Figure 4 shows the time dependence of the L ²-error

$\Vert |{\psi }_N\left(t\right)\rangle -|{\psi }_{\text{ex}}\left(t\right)\rangle {\Vert }_{L^2},(59)$

where the semiclassical wavefunctions were generated using N = 2¹⁶ trajectories. Sampling from the square root of the Husimi density (Case sqrt-H) results in approximately twice smaller error than sampling from the Husimi density itself. For (Case sqrt-H), the figure also displays the analytical error estimate derived in (55), which matches the numerical error up to small statistical noise. To remove this statistical noise and to match the analytical estimate (55) more accurately, in the lower panel of Figure 4 we plot the empirical root mean square error (RMSE) [59] S ₁₀₀, where

$S_K:=\frac{1}{K}\sum _{j=1}^K\Vert |{\psi }_N^{\left(j\right)}\left(t\right)\rangle -|{\psi }_{\text{ex}}\left(t\right)\rangle {\Vert }^2,(60)$

K is an integer number of independent simulations (indexed by j) and each ${\psi }_N^{(j)}$ was itself approximated using N independent samples. One can see that

$\mathbb{E}\left[S_K\right]=\mathbb{E}\left[\Vert |{\psi }_N\left(t\right)\rangle -|{\psi }_{\text{ex}}\left(t\right)\rangle {\Vert }^2\right]=\frac{\mathbb{V}\left[|\tilde{\psi }\left(t\right)\rangle \right]}{N}.(61)$

For K = 1, one obtains the result represented in the upper panel of Figure 4. For K → ∞, due to the strong law of large numbers, S _K converges almost surely to its expectation value and hence to the analytical error estimate.

FIGURE 4

FIGURE 4. Time dependence of the sampling error of the Herman–Kluk wavefunction propagated in a harmonic oscillator. The upper panel is produced by one run with N = 2¹⁶ = 65,536 trajectories, whereas the lower panel is produced by K = 100 independent runs, each with N = 2¹⁶ = 65,536 trajectories, and averaging the square of the error over the K runs. The analytical error estimate for the sampling from the square root of the Husimi density (Case sqrt-H) is shown with the dotted line.

The lower panel of Figure 4 shows the empirical RMSE computed with K = 100 and N = 2¹⁶. For (Case sqrt-H), it coincides almost perfectly with the analytical prediction. We note that the error reaches its maximum whenever the wavefunction passes through the bottom of the potential and its minimum when the wavefunction arrives at the turning points. Even though this analysis makes only sense for finite variance, the figure also displays the empirical RMSE for (Case H), which is nearly constant. Both panels show a qualitatively similar error evolution for (Case H), however, at a magnitude that is considerably larger than the error for (Case sqrt-H).

4.3 Morse potential

To investigate the convergence of the Herman–Kluk wavefunction in anharmonic systems, we consider dynamics generated by a less and more anharmonic Morse potentials. The parameters were taken from [29]. Our initial state is a Gaussian wavepacket with zero initial position and momentum (q ₀, p ₀) = (0, 0), and with a width parameter γ = 0.00456 a.u. $\approx 1000{\text{cm}}^{-1}$ . The Morse potential

$V\left(x\right)=V_{\text{eq}}+D_e{\left[1-e^{-a\left(x-q_{\text{eq}}\right)}\right]}^2(62)$

is characterized by the dissociation energy D _e , decay parameter a, and the position q _eq and energy V _eq of the minimum. We considered two Morse potentials, both with V _eq = 0.1 and q _eq = 20.95 a.u., but with different values of a and D _e . The latter two parameters, however, were chosen so that the global harmonic potential fitted to the Morse potential at q _eq had the same frequency

${\omega }_{\text{eq}}=\sqrt{V^{\prime \prime }\left(q_{\text{eq}}\right)}=\sqrt{2D_e{a}^2}=0.0041\text{\hspace{0.17em}a.u.}\approx 900{\text{cm}}^{-1}(63)$

for the two Morse potentials. The anharmonicity of the potentials was conveniently controlled with the dimensionless parameter

$\chi =\hslash \frac{{\omega }_{\text{eq}}}{4D_e},(64)$

which is also reflected in the bound energy levels [71]

$E_n=\hslash {\omega }_{\text{eq}}\left[\left(n+\frac{1}{2}\right)-\chi {\left(n+\frac{1}{2}\right)}^2\right](65)$

of a Morse oscillator. Then D _e and a are given by

$\begin{array}{cc}\hfill D_e& =\hslash \frac{{\omega }_{\text{eq}}}{4\chi },\hfill \\ \hfill a& =\sqrt{2\frac{{\omega }_{\text{eq}}\chi }{\hslash }}.\hfill \end{array}(66)$

We choose two different values

$\chi =0.005\text{and}\chi =0.01(67)$

of anharmonicity and compare the Herman–Kluk propagation with a grid-based reference quantum calculation obtained by the Fourier-split method [72], which is second-order accurate with respect to the time step. The position grid was set from x = −200 to 1500 with 4096 equidistant points for χ = 0.005 and from x = −200 to 10,000 with 16,384 equidistant points for χ = 0.01. The larger grid for χ = 0.01 was required to capture oscillations of the wavefunctions in the tail region, which are due to the increased anharmonicity. For the time propagation of the Herman–Kluk wavefunction, we used a second order Størmer-Verlet scheme with step size Δt = 8 a.u. $\approx 0.194$ fs.

Because the Herman–Kluk approximation is not exact in a Morse potential, to separate the statistical convergence error from the semiclassical error of the fully converged Herman–Kluk approximation, in Figure 5 we show the L ²-error

$\Vert {\psi }_N\left(t\right)-{\psi }_{2N}\left(t\right)\Vert (68)$

between the Herman–Kluk wavefunctions calculated with N and 2N trajectories as a function of N. Both wavefunctions were propagated in a Morse potential with anharmonicity χ = 0.005. Each panel includes the error for the fixed time t after approximately one oscillation (196 time steps) and ten oscillations (1960 time steps). In addition, the convergence rates for both sampling schemes were fitted to the same power law (57). We observe that sampling from the square root of the Husimi density (Case sqrt-H), shown in the lower panel, performs better than sampling from the Husimi density (Case H), displayed in the upper panel.

FIGURE 5

FIGURE 5. Sampling error between the Herman–Kluk wavefunctions obtained with N and 2N Monte Carlo quadrature points as a function of N. The wavefunctions are calculated in a Morse potential with anharmonicity parameter χ = 0.005 after approximately one (solid line) or ten oscillations (dashed line). The upper panel shows both errors and their approximated convergence rates for (Case H). Similarly, (Case sqrt-H) is displayed in the lower panel.

Figure 6 shows the analogous results obtained in a Morse potential with a larger anharmonicity parameter χ = 0.01. Here, we display the wavefunctions after 202 and 2020 time steps, which are approximately the times after the first and tenth oscillations. As expected for anharmonic evolution, the error after ten oscillations is worse than after one oscillation. Increased anharmonicity also increases the error in comparison to Figure 5. Again, the sampling from the square root of the Husimi density (Case sqrt-H) has consistently a lower error than (Case H).

FIGURE 6

FIGURE 6. Same as Figure 5, except that the anharmonicity parameter of the Morse potential was increased to χ = 0.01.

To complement the abstract convergence study of the L ²-error, in the upper panels of Figures 7 and 8, we compare the more intuitive position probability densities of the “exact” quantum grid-based solution with those of (Case H) and (Case sqrt-H). Both figures were obtained in a Morse potential with anharmonicity parameter χ = 0.01 at a time after ten oscillations. The difference lies in the number of trajectories used. The less converged results in Figure 7 were obtained with only N = 800 trajectories, whereas the fully converged results in Figure 8 employed N = 100 ⋅ 2¹⁴ trajectories. The lower panels of the two figures display the absolute errors of the position density for the two cases, measured with respect to the “exact” grid-based position density. The two panels of Figure 7 confirm again that sampling from the square root of Husimi density (Case sqrt-H) results in faster convergence than sampling from the Husimi density (Case sqrt-H). The upper panel of Figure 8 shows that in the numerically converged regime, the Herman–Kluk propagator approximates the exact solution in this system very well, regardless whether one samples from the Husimi density or its square root. The fact that results are numerically converged is confirmed in the lower panel of Figure 8, where the errors of (Case H) and (Case sqrt-H) are approximately the same, which implies that the common remaining error is the error of the Herman–Kluk approximation and not the phase-space discretization error.

FIGURE 7

FIGURE 7. Position probability densities (upper panel) and their absolute errors (lower panel) in a Morse potential with χ = 0.01 after 10 oscillations (2020 time steps $\approx 392$ fs). Position probability densities of the Herman–Kluk approximation for (Case H) and (Case sqrt-H) were computed with N = 800 trajectories.

FIGURE 8

FIGURE 8. Same as Figure 7, except that the results were computed with N = 100 ⋅ 2¹⁴ ≈ 1.6 ⋅ 10⁶ trajectories and, therefore, are converged.

Analytical expressions and numerical fits to the convergence rates for all studied systems are summarized in Table 1.

TABLE 1

TABLE 1. Summary of convergence rates for (Case sqrt-H) and (Case H) at initial time, in a harmonic oscillator and in two Morse potentials after one and ten oscillations.

Finally, we note that Figure 1 in the Introduction was based, as Figure 8, on the Morse potential with anharmonicity parameter χ = 0.01 and Herman–Kluk calculations used N = 100 ⋅ 2¹⁴ trajectories. For the computation of the position density, we used the more efficient (Case sqrt-H). The spectra in the upper panel were obtained by the Fourier transform of the wavepacket autocorrelation function; the Herman–Kluk autocorrelation function was computed by sampling from the Husimi density (Case H), because it gives the exact result (=1) at t = 0 regardless of the number of trajectories and, therefore, converges faster at short times.

5 Conclusion and outlook

We compared two different sampling strategies for evaluating the semiclassical wavefunction evolved with the Herman–Kluk propagator. For the initial phase-space sampling, we either used the Husimi density (Case H) or its square root (Case sqrt-H). We showed that the square root sampling produces a Monte Carlo integrand with finite second moment, while the Husimi sampling comes with an undesirable infinite second moment. The numerical experiments for the harmonic oscillator and two Morse oscillators with different extents of anharmonicity confirm that the infinite second moment results in a slower convergence of the Monte Carlo estimator. Therefore, we recommend the square root approach (Case sqrt-H) whenever the Herman–Kluk propagator is used directly for approximating the quantum-mechanical wavefunction and the L ²-error of the wavefunction approximation is the relevant accuracy measure. However, we explicitly verified that at initial time the square root density’s second moment, even though it is finite, has an unfavorable, exponential dependence on the dimension, possibly leading to a large number of trajectories required for a reasonably accurate wavefunction.

Although the wavefunction is a central object in quantum mechanics, one is often interested directly in observables, which can be computed as expectation values. It is clearly inefficient to compute expectation values, such as energy or squared norm, with the Herman–Kluk propagator by computing the wave function first. A follow-up paper, in which an analysis similar to that presented here will be applied to the autocorrelation function as well as to the expectation values is in preparation. In particular, the approach to expectation values proposed in [52] will be analysed in detail. Our present analysis and sampling approaches could, in principle, help increase the efficiency of any Gaussian-based method, although it is difficult to predict the possible implications that the coupling between different Gaussians present in Gaussian basis methods [9, 10] might have for the choice of the optimal sampling density.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Author contributions

All authors have contributed to the design and execution of the research and to the writing of the manuscript. FK performed all numerical calculations and prepared all figures.

Acknowledgments

The authors acknowledge financial support from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (Grant Agreement No. 683069–MOLEQULE) and from the EPFL.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphy.2023.1106324/full#supplementary-material

References

1. Miller WH. Classical S matrix: Numerical application to inelastic collisions. J Chem Phys (1970) 53:3578–87. doi:10.1063/1.1674535

CrossRef Full Text | Google Scholar

2. Miller WH. The semiclassical initial value representation: A potentially practical way for adding quantum effects to classical molecular dynamics simulations. J Phys Chem A (2001) 105:2942–55. doi:10.1021/jp003712k

CrossRef Full Text | Google Scholar

3. Heller EJ. The semiclassical way to dynamics and spectroscopy. Princeton, NJ: Princeton University Press (2018).

Google Scholar

4. Heller EJ. Frozen Gaussians: A very simple semiclassical approximation. J Chem Phys (1981) 75:2923–31. doi:10.1063/1.442382

CrossRef Full Text | Google Scholar

5. Heller EJ. Time-dependent approach to semiclassical dynamics. J Chem Phys (1975) 62:1544–55. doi:10.1063/1.430620

CrossRef Full Text | Google Scholar

6. Herman MF, Kluk E. A semiclasical justification for the use of non-spreading wavepackets in dynamics calculations. Chem Phys (1984) 91:27–34. doi:10.1016/0301-0104(84)80039-7

CrossRef Full Text | Google Scholar

7. Herman MF. Time reversal and unitarity in the frozen Gaussian approximation for semiclassical scattering. J Chem Phys (1986) 85:2069–76. doi:10.1063/1.451150

CrossRef Full Text | Google Scholar

8. Kluk E, Herman MF, Davis HL. Comparison of the propagation of semiclassical frozen Gaussian wave functions with quantum propagation for a highly excited anharmonic oscillator. J Chem Phys (1986) 84:326–34. doi:10.1063/1.450142

CrossRef Full Text | Google Scholar

9. Shalashilin DV, Child MS. The phase space ccs approach to quantum and semiclassical molecular dynamics for high-dimensional systems. Chem Phys (2004) 304:103–20. doi:10.1016/j.chemphys.2004.06.013

CrossRef Full Text | Google Scholar

10. Richings G, Polyak I, Spinlove K, Worth G, Burghardt I, Lasorne B. Quantum dynamics simulations using Gaussian wavepackets: The vmcg method. Int Rev Phys Chem (2015) 34:269–308. doi:10.1080/0144235X.2015.1051354

CrossRef Full Text | Google Scholar

11. Kay KG. Semiclassical initial value treatments of atoms and molecules. Annu Rev Phys Chem (2005) 56:255–80. doi:10.1146/annurev.physchem.56.092503.141257

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Kay KG. Numerical study of semiclassical initial value methods for dynamics. J Chem Phys (1994) 100:4432–45. doi:10.1063/1.466273

CrossRef Full Text | Google Scholar

13. Walton AR, Manolopoulos DE. A new semiclassical initial value method for franck-condon spectra. Mol Phys (1996) 87:961–78. doi:10.1080/00268979600100651

CrossRef Full Text | Google Scholar

14. Garashchuk S, Tannor D. Wave packet correlation function approach to H2(ν)+H→H+H2(ν′): Semiclassical implementation. Chem Phys Lett (1996) 262:477–85. doi:10.1016/0009-2614(96)01111-6

CrossRef Full Text | Google Scholar

15. Thoss M, Wang H. Semiclassical description of molecular dynamics based on initial-value representation methods. Annu Rev Phys Chem (2004) 55:299–332. doi:10.1146/annurev.physchem.55.091602.094429

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Spanner M, Batista VS, Brumer P. Is the filinov integral conditioning technique useful in semiclassical initial value representation methods? J Chem Phys (2005) 122:084111. doi:10.1063/1.1854634

CrossRef Full Text | Google Scholar

17. Tatchen J, Pollak E. Semiclassical on-the-fly computation of the S ₀ → S ₁ absorption spectrum of formaldehyde. J Chem Phys (2009) 130:041103. doi:10.1063/1.3074100

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Ceotto M, Atahan S, Shim S, Tantardini GF, Aspuru-Guzik A. First-principles semiclassical initial value representation molecular dynamics. Phys Chem Chem Phys (2009) 11:3861–7. doi:10.1039/B820785B

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Kay KG. Integral expressions for the semiclassical time-dependent propagator. J Chem Phys (1994) 100:4377–92. doi:10.1063/1.466320

CrossRef Full Text | Google Scholar

20. Miller WH. On the relation between the semiclassical initial value representation and an exact quantum expansion in time-dependent coherent states. J Phys Chem B (2002) 106:8132–5. doi:10.1021/jp020500+

CrossRef Full Text | Google Scholar

21. Miller WH. An alternate derivation of the Herman–Kluk (coherent state) semiclassical initial value representation of the time evolution operator. Mol Phys (2002) 100:397–400. doi:10.1080/00268970110069029

CrossRef Full Text | Google Scholar

22. Deshpande SA, Ezra GS. On the derivation of the herman–kluk propagator. J Phys A (2006) 39:5067–78. doi:10.1088/0305-4470/39/18/020

CrossRef Full Text | Google Scholar

23. Tannor DJ, Heller EJ. Polyatomic Raman scattering for general harmonic potentials. J Chem Phys (1982) 77:202–18. doi:10.1063/1.443643

CrossRef Full Text | Google Scholar

24. Begušić T, Roulet J, Vaníček J. On-the-fly ab initio semiclassical evaluation of time-resolved electronic spectra. J Chem Phys (2018) 149:244115. doi:10.1063/1.5054586

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Hagedorn GA. Semiclassical quantum mechanics. I. The ℏ→ 0 limit for coherent states. Commun Math Phys (1980) 71:77–93. doi:10.1007/bf01230088

CrossRef Full Text | Google Scholar

26. Lee SY, Heller EJ. Exact time-dependent wave packet propagation: Application to the photodissociation of methyl iodide. J Chem Phys (1982) 76:3035–44. doi:10.1063/1.443342

CrossRef Full Text | Google Scholar

27. Hagedorn GA. Raising and lowering operators for semiclassical wave packets. Ann Phys (Ny) (1998) 269:77–104. doi:10.1006/aphy.1998.5843

CrossRef Full Text | Google Scholar

28. Coalson RD, Karplus M. Multidimensional variational Gaussian wave packet dynamics with application to photodissociation spectroscopy. J Chem Phys (1990) 93:3919–30. doi:10.1063/1.458778

CrossRef Full Text | Google Scholar

29. Begušić T, Cordova M, Vaníček J. Single-Hessian thawed Gaussian approximation. J Chem Phys (2019) 150:154117. doi:10.1063/1.5090122

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Prlj A, Begušić T, Zhang ZT, Fish GC, Wehrle M, Zimmermann T, et al. Semiclassical approach to photophysics beyond kasha’s rule and vibronic spectroscopy beyond the condon approximation. The case of azulene. J Chem Theor Comput. (2020) 16:2617–26. doi:10.1021/acs.jctc.0c00079

CrossRef Full Text | Google Scholar

31. Wu Y, Herman MF. Nonadiabatic surface hopping herman-kluk semiclassical initial value representation method revisited: Applications to tully’s three model systems. J Chem Phys (2005) 123:144106. doi:10.1063/1.2049251

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Wu Y, Herman MF. A justification for a nonadiabatic surface hopping herman-kluk semiclassical initial value representation of the time evolution operator. J Chem Phys (2006) 125:154116. doi:10.1063/1.2358352

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Hochman G, Kay KG. Tunneling by a semiclassical initial value method with higher order corrections. J Phys A (2008) 41:385303. doi:10.1088/1751-8113/41/38/385303

CrossRef Full Text | Google Scholar

34. Hochman G, Kay KG. Tunneling in two-dimensional systems using a higher-order herman–kluk approximation. J Chem Phys (2009) 130:061104. doi:10.1063/1.3079544

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Filinov VS. Calculation of the feynman integrals by means of the Monte Carlo method. Nucl Phys B (1986) 271:717–25. doi:10.1016/S0550-3213(86)80034-7

CrossRef Full Text | Google Scholar

36. Makri N, Miller WH. Monte Carlo integration with oscillatory integrands: Implications for feynman path integration in real time. Chem Phys Lett (1987) 139:10–4. doi:10.1016/0009-2614(87)80142-2

CrossRef Full Text | Google Scholar

37. Makri N, Miller WH. Monte Carlo path integration for the real time propagator. J Chem Phys (1988) 89:2170–7. doi:10.1063/1.455061

CrossRef Full Text | Google Scholar

38. Elran Y, Kay KG. Improving the efficiency of the herman–kluk propagator by time integration. J Chem Phys (1999) 110:3653–9. doi:10.1063/1.478255

CrossRef Full Text | Google Scholar

39. Kaledin AL, Miller WH. Time averaging the semiclassical initial value representation for the calculation of vibrational energy levels. J Chem Phys (2003) 118:7174–82. doi:10.1063/1.1562158

CrossRef Full Text | Google Scholar

40. Buchholz M, Grossmann F, Ceotto M. Mixed semiclassical initial value representation time-averaging propagator for spectroscopic calculations. J Chem Phys (2016) 144:094102. doi:10.1063/1.4942536

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Buchholz M, Grossmann F, Ceotto M. Simplified approach to the mixed time-averaging semiclassical initial value representation for the calculation of dense vibrational spectra. J Chem Phys (2018) 148:114107. doi:10.1063/1.5020144

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Shao J, Makri N. Forward-backward semiclassical dynamics in the interaction representation. J Chem Phys (2000) 113:3681–5. doi:10.1063/1.1287823

CrossRef Full Text | Google Scholar

43. Petersen J, Pollak E. Semiclassical initial value representation for the quantum propagator in the heisenberg interaction representation. J Chem Phys (2015) 143:224114. doi:10.1063/1.4936922

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Ceotto M, Atahan S, Tantardini GF, Aspuru-Guzik A. Multiple coherent states for first-principles semiclassical initial value representation molecular dynamics. J Chem Phys (2009) 130:234113. doi:10.1063/1.3155062

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Grossmann F. A semiclassical hybrid approach to many particle quantum dynamics. J Chem Phys (2006) 125:014111. doi:10.1063/1.2213255

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Goletz CM, Grossmann F. Decoherence and dissipation in a molecular system coupled to an environment: An application of semiclassical hybrid dynamics. J Chem Phys (2009) 130:244107. doi:10.1063/1.3157162

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Grossmann F. A semiclassical hybrid approach to linear response functions for infrared spectroscopy. Phys Scr T (2016) 91:044004. doi:10.1088/0031-8949/91/4/044004

CrossRef Full Text | Google Scholar

48. Antipov SV, Ye Z, Ananth N. Dynamically consistent method for mixed quantum-classical simulations: A semiclassical approach. J Chem Phys (2015) 142:184102. doi:10.1063/1.4919667

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Church MS, Antipov SV, Ananth N. Validating and implementing modified filinov phase filtration in semiclassical dynamics. J Chem Phys (2017) 146:234104. doi:10.1063/1.4986645

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Malpathak S, Church MS, Ananth N. A semiclassical framework for mixed quantum classical dynamics. J Phys Chem A (2022) 126:6359–75. doi:10.1021/acs.jpca.2c03467

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Pollak E, Upadhyayula S, Liu J. Coherent state representation of thermal correlation functions with applications to rate theory. J Chem Phys (2022) 156:244101. doi:10.1063/5.0088163

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Lasser C, Sattlegger D. Discretising the herman-kluk propagator. Numerische Mathematik (2017) 137:119–57. doi:10.1007/s00211-017-0871-0

CrossRef Full Text | Google Scholar

53. Lasser C, Lubich C. Computing quantum dynamics in the semiclassical regime. Acta Numerica (2020) 29:229–401. doi:10.1017/S0962492920000033

CrossRef Full Text | Google Scholar

54. Hall BC. Quantum theory for mathematicians. New York: Springer (2013). doi:10.1007/978-1-4614-7116-5

CrossRef Full Text | Google Scholar

55. Martinez A. An introduction to semiclassical and microlocal analysis. New York: Springer-Verlag (2002).

Google Scholar

56. Kay KG. The herman–kluk approximation: Derivation and semiclassical corrections. Chem Phys (2006) 322:3–12. doi:10.1016/j.chemphys.2005.06.019

CrossRef Full Text | Google Scholar

57. Swart T, Rousse V. A mathematical justification for the herman-kluk propagator. Commun Math Phys (2008) 286:725–50. doi:10.1007/s00220-008-0681-4

CrossRef Full Text | Google Scholar

58. Robert D. On the herman–kluk semiclassical approximation. Rev Math Phys (2010) 22:1123–45. doi:10.1142/s0129055x1000417x

CrossRef Full Text | Google Scholar

59. Caflisch RE. Monte Carlo and quasi-Monte Carlo methods. Acta Numerica (1998) 7:1–49. doi:10.1017/S0962492900002804

CrossRef Full Text | Google Scholar

60. Brewer ML, Hulme JS, Manolopoulos DE. Semiclassical dynamics in up to 15 coupled vibrational degrees of freedom. J Chem Phys (1997) 106:4832–9. doi:10.1063/1.473532

CrossRef Full Text | Google Scholar

61. Hairer E, Lubich C, Wanner G. Geometric numerical integration: Structure-preserving algorithms for ordinary differential equations. New York: Springer Berlin Heidelberg (2006).

Google Scholar

62. Hairer E, Lubich C, Wanner G. Geometric numerical integration illustrated by the Størmer-Verlet method. Acta Numerica (2003) 12:399–450. doi:10.1017/S0962492902000144

CrossRef Full Text | Google Scholar

63. Stewart I, Tall D. Complex analysis. Cambridge, United Kingdom: Cambridge University Press (2018). doi:10.1017/9781108505468

CrossRef Full Text | Google Scholar

64. Kahan W, Li RC. Composition constants for raising the orders of unconventional schemes for ordinary differential equations. Math Comput (1997) 66:1089–99. doi:10.1090/s0025-5718-97-00873-9

CrossRef Full Text | Google Scholar

65. Suzuki M. General theory of higher-order decomposition of exponential operators and symplectic integrators. Phys Lett A (1992) 165:387–95. doi:10.1016/0375-9601(92)90335-J

CrossRef Full Text | Google Scholar

66. Vaníček J. Dephasing representation of quantum fidelity for general pure and mixed states. Phys Rev E (2006) 73:046204. doi:10.1103/PhysRevE.73.046204

CrossRef Full Text | Google Scholar

67. Mollica C, Vaníček J. Beating the efficiency of both quantum and classical simulations with a semiclassical method. Phys Rev Lett (2011) 107:214101. doi:10.1103/PhysRevLett.107.214101

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Mollica C, Zimmermann T, Vaníček J. Efficient sampling avoids the exponential wall in classical simulations of fidelity. Phys Rev E (2011) 84:066205. doi:10.1103/PhysRevE.84.066205

CrossRef Full Text | Google Scholar

69. Zimmermann T, Vaníček J. Role of sampling in evaluating classical time autocorrelation functions. J Chem Phys (2013) 139:104105. doi:10.1063/1.4820420

PubMed Abstract | CrossRef Full Text | Google Scholar

70. Durrett R. Probability: Theory and examples. Cambridge, United Kingdom: Cambridge University Press (2019). doi:10.1017/9781108591034

CrossRef Full Text | Google Scholar

71. Tannor DJ. Introduction to quantum mechanics: A time-dependent perspective. Sausalito: University Science Books (2007).

Google Scholar

72. Feit MD, Fleck JA. Solution of the Schrödinger equation by a spectral method ii: Vibrational energy levels of triatomic molecules. J Chem Phys (1983) 78:301–8. doi:10.1063/1.444501

CrossRef Full Text | Google Scholar

73. Dyke P. An introduction to laplace transforms and fourier series. London: Springer London (2013). doi:10.1007/978-1-4471-6395-4

CrossRef Full Text | Google Scholar