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.


I. INTRODUCTION
Semiclassical initial value representation techniques 1,2 have evolved into useful tools for calculations of the dynamics of atoms and molecules. 3Frozen Gaussians and their superposition were introduced 4 by Heller in 1981 as an extension to the thawed Gaussian approximation 5 in order to capture nonlinear spreading of the wavepackets.Herman and Kluk justified the frozen-Gaussian ansatz and introduced a refined approximation, 6,7 now known as the Herman-Kluk propagator, which contains an additional prefactor that rigorously compensates for the fixed width of these frozen Gaussians.
In the spirit of the semiclassical initial value representation, the Herman-Kluk propagator avoids the root search problem. 2,84][25][26][27][28] In chaotic and other systems, however, increased anharmonicity leads to wavepacket splitting and nontrivial 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., Fig. 1).
The multi-trajectory nature of the Herman-Kluk propagator is not only its advantage, but also its bottleneck.
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 methods, the advantage of this approach is that the converged results agree exactly with the converged results of the original Herman-Kluk propagator.For this reason, 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 of the initial state. 2,45In contrast, 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.
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 Sec. 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 Sec. 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.

II. DISCRETISING THE HERMAN-KLUK PROPAGATOR A. Herman-Kluk propagator
Evolution of a quantum state |ψ(t) is governed by the time-dependent Schrödinger equation where is the reduced Planck constant and Ĥ is the Hamiltonian operator.Here, we assume the Hamiltonian to be the time-independent operator 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 spectraltheorem 46 provides for all times t ∈ R a well-defined unitary propagator in terms of which the solution of (1) can be expressed as for all square integrable initial data 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 −1 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 with a centre at the phase-space point z = (q, p) ∈ R 2D and with a fixed, real, symmetric, positive-definite widthmatrix γ ∈ R D×D .The frozen Gaussians {|g γ z : z ∈ R 2D } form an over-complete subset of the Hilbert space of square integrable functions.They have the striking property that any |ψ ∈ L 2 (R D ) can be decomposed as with the scaled phase-space measure dν = dz/(2π ) D .Using this decomposition for our solution of the timedependent Schrödinger equation, we obtain Based on the continuous superposition (7), the state obtained by applying the Herman-Kluk propagator to |ψ 0 can be expressed as where z(t) = (q(t), p(t)) is the solution to the underlying classical Hamiltonian system for the Hamiltonian function h(q, p) = T (p)+V (q), where is the standard symplectic matrix.S denotes the classical action integral which solves the initial value problem for all z = (q, p) ∈ R 2D .In ( 8), the Herman-Kluk prefactor depends on the matrices 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 with Hess h the Hessian of the Hamiltonian function h.9][50][51] It has been shown that the exactly evolved quantum state is approximated by the state obtained by applying the Herman-Kluk propagator to the initial state |ψ 0 , with an error of the order of .More precisely, for all initial data |ψ 0 of norm one, where T > 0 is a fixed time and C(T ) > 0 is a constant independent of .If the potential is at most harmonic, then the approximation is exact. 49

B. Discretisation
Evolving the wavefunction (8) with the Herman-Kluk propagator requires evaluating an integral over the phase space 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 ∈ 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 R 2D must be evaluated using mesh-free discretization, such as Monte Carlo methods. 52o evaluate the Herman-Kluk wavefunction by Monte Carlo sampling, we rewrite the integral from Eq. ( 8) as where we introduced the notation 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 i.e., a propagated frozen Gaussian multiplied by the Herman-Kluk and phase prefactors, and a timeindependent function The Monte Carlo estimator is then given by where ψ(t, z j ) is the state obtained from initial condition z j and z 1 , z 2 , . . ., 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. ( 8) into the product g γ z |ψ 0 = r(z)ρ(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 of the initial state |ψ 0 .For Husimi distribution, the prefactor r(z) becomes However, since the wavefunction evolved with the Herman-Kluk propagator is first-and not second-order in |ψ 0 , it is natural to also test the square root of the Husimi density and to consider the probability density with a bounded prefactor The two probability densities ρ H (z) and ρ(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 In general, the integral over 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 |ψ 0 = |g γ z0 centred at some phase-space point z 0 = (q 0 , p 0 ) ∈ R 2D , then where Σ 0 := γ 0 0 γ −1 is the matrix containing the width parameters of the initial coherent state.Then, the Husimi density is given by whereas the second approach provides the density and prefactor C. Summary of the numerical algorithm Taking into account all the previous considerations, we slightly extend the natural numerical algorithm (described, e.g., in Ref. 53, Sec.4]) for finding the Herman-Kluk approximation to the wavefunction at time t: Algorithm 1. (Herman-Kluk propagation) 1. Draw independent samples z 1 , . . ., z N ∈ R 2D from a distribution with density ρ H (z) or ρ(z) given by Eqs. ( 28) and (29).
We note that Eqs. ( 9), ( 12) and ( 15) can be evaluated simultaneously with a single numerical integrator.To increase the accuracy of the time integration one can use higher-order composition methods. 58,59They 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.

III. PHASE SPACE ERROR ANALYSIS
Algorithm 1 relies on the discretisation of the phasespace integrals and of the system of ordinary differential equations.Here we focus only on the phase-space discretisation error.

A. Moments of the integrand
To assess the accuracy of the Monte Carlo estimator (22), we examine the moments of the two different integrands: |ψ H (t) and | ψ(t) .First, we observe that for any s > 0, (33)   since 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 Appendix A. In contrast, 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. 60) provides convergence of the estimator, with probability one.The non-existing second moment for (Case H) results in a slightly worse convergence rate than the one for (Case sqrt-H).Numerical results in Sec.
IV confirm this expectation.

B. Mean squared error
For (Case sqrt-H) the second moment is finite, so that the mean squared error of the Monte Carlo estimator (22)  is well-defined and satisfies where the expectation value and the variance are with respect to the density ρ(z).Moreover, 53 In the special case of an initial Gaussian initial wavepacket |ψ 0 = |g γ z0 , this simplifies to and at initial time t = 0 we obtain since the Herman-Kluk prefactor satisfies R(0, z) = 1.
For a more general assessment of the variance, numerical experiments in Sec.IV 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 Note that Eq. ( 36) 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 Appendix B, we explain how Eq. ( 36) can also provide a rigorous lower bound and asymptotic estimate of the number of trajectories needed for convergence of a single simulation.

C. Other sampling densities
For the special case of an initial Gaussian wavepacket |ψ 0 = |g γ z0 , the two proposed sampling densities ρ H (z) and ρ(z) belong to a family of normal distributions with probability density functions with a ≥ 2. In the spirit of importance sampling, the Herman-Kluk wavefunction can accordingly be written as a phase-space average At time t = 0, the norm of the integrand satisfies which implies for the variance For a > 2, the function a → a 2 /(2a − 4) attains its minimum at a = 4, which corresponds to the sampling density ρ.In other words, (Case sqrt-H) is optimal as far as the mean squared error is concerned.

D. Harmonic motion
The one-dimensional harmonic oscillator 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 |ψ 0 = |g γ z0 with position and momentum q 0 , p 0 ∈ R, the exact wavefunction is given by 61 where p t = p 0 cos(ωt) − bq 0 sin(ωt) and ( 49) with the abbreviations z t = b cos(ωt) + α 0 sin(ωt), α 0 = iγ and b = mω.Here, β 0 includes the normalization constant for the wavefunction at time t = 0.The classical action is given by The four components of stability matrix M can be obtained, from their definition (14), by differentiating expressions for q t and p t with respect to q 0 and p 0 , namely The Herman-Kluk prefactor satisfies so that the variance (37) can be written as 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

IV. 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 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. ( 56) 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 phasespace centre z 0 ∈ R 2D .

A. Initial phase-space error
We start by considering a spherical initial Gaussian wavepacket |ψ ex (0) = |g γ z0 with a width parameter γ = 2Id D in one and four dimensions.For D = 1, it is centreed at z 0 = (−1, 0) and for D = 4, at z 0 = −(1, 1, 1, 1, 0, 0, 0, 0). Figure 2 shows the L 2 -distance 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 (39) of the mean squared error ( 36) of (Case sqrt-H) is fulfilled.For (Case H), the curve-fitting approximation ( 56

B. Harmonic potential
To analyze the effect of dynamics on the convergence, let us consider a harmonic potential V (x) = x 2 /2 in one dimension and explore the dynamics for one full oscillation period.The initial Gaussian wavepacket is still localized in z 0 = (−1, 0) with width γ = 2.In Fig. 4, we compare the exact wavefunction (46) 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. ( 48), ( 49), ( 51) and (53).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 Fig. 4 shows the time dependence of the L 2 -error where the semiclassical wavefunctions were generated using N = 2 16 trajectories.Sampling from the square root of the Husimi density (Case sqrt-H) results in ap-proximately twice smaller error than sampling from the Husimi density itself.For (Case sqrt-H), the figure also displays the analytical error estimate derived in (54), which matches the numerical error up to small statistical noise.To remove this statistical noise and to match the analytical estimate (54) more accurately, in the lower panel of Fig. 4 we plot the empirical root mean square error (RMSE) 52 S 100 , where K is an integer number of independent simulations (indexed by j) and each ψ N was itself approximated using N independent samples.One can see that For K = 1, one obtains the result represented in the upper panel of Fig. 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.The lower panel of Fig. 4 shows the empirical RMSE computed with K = 100 and N = 2 16 .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).

C. 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 [27].Our initial state is a Gaussian wavepacket with zero initial position and momentum, (q 0 , p 0 ) = (0, 0), and with a width parameter γ = 0.00456 a.u.≈ 1000 cm −1 .The Morse potential 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.95a.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 ω eq = V ′′ (q eq ) = 2D e a 2 = 0.0041 a.u.≈ 900 cm −1 (62)   for the two Morse potentials.The anharmonicity of the potentials was conveniently controlled with the dimensionless parameter which is also reflected in the bound energy levels 61  timesteps) and ten oscillations (1960 timesteps).In addition, the convergence rates for both sampling schemes were fitted to the same power law (56).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 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 Fig. 5. Again, the sampling from the square root of the Husimi density (Case sqrt-H) has consistently a lower error than (Case H).
To complement the abstract convergence study of the L 2 -error, in the upper panels of Fig. 7 and Fig. 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 param- eter χ = 0.01 at a time after ten oscillations.The difference lies in the number of trajectories used.The less converged results in Fig. 7 were obtained with only N = 800 trajectories, whereas the fully converged results in Fig. 8 employed N = 100 • 2 14 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 Fig. 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 Fig. 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 Fig. 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.
Analytical expressions and numerical fits to the convergence rates for all studied systems are summarized in Table I.
Finally, we note that Fig. 1    monicity parameter χ = 0.01 and Herman-Kluk calculations used N = 100 • 2 14 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.

V. 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 2 -error of the wavefunction approximation is the relevant accuracy measure.A followup paper applying a similar analysis to the autocorrelation function as well as to the expectation values computed with the Herman-Kluk propagator is in preparation.

FIG. 1 .
FIG.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 Sec.IV.
) 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, Fig.3displays the L 2 -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 5 trajectories.The analytical prediction(39) for (Case sqrt-H) is realized.Moreover, the error for (Case H) increases faster with D.

4 D 36 FIG. 2 .
FIG.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).

4 D − 1 • N −1/ 2 FIG. 3 .
FIG. 3. Dependence of the sampling error of the initial wavefunction on dimension D for N = 100 • 2 13 ≈ 8 • 10 5 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.

FIG. 4 .
FIG.4.Time dependence of the sampling error of the Herman-Kluk wavefunction propagated in a harmonic oscillator.The upper panel is produced by one independent run with N = 216 = 65536 trajectories, whereas the lower panel is produced by K = 100 independent runs, each with N = 216 = 65536 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.

2 HusimiFIG. 5 .
FIG. 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.
Morse oscillator.Then D e and a are given by D e = ω eq 4χ , a = 2 ω eq χ .

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