Validating and Implementing Modified Filinov Phase Filtration in Semiclassical Dynamics

TL;DR

This paper demonstrates that the modified Filinov filtration significantly improves the efficiency and accuracy of semiclassical correlation function calculations in MQC-IVR, especially by damping noise in oscillatory regions and enabling linear scaling.

Contribution

It introduces a novel MQC-IVR formulation with linear scaling and validates the effectiveness of MFF in reducing noise and computational effort in semiclassical dynamics.

Findings

MFF reduces noise by damping amplitude in oscillatory phase regions.

The new MQC-IVR formulation achieves linear computational scaling.

MFF is highly effective for correlation function calculations, unlike in wavepacket propagation.

Abstract

The Mixed Quantum-Classical Initial Value Representation (MQC-IVR) is a recently introduced approximate semiclassical (SC) method for the calculation of real-time quantum correlation functions. MQC-IVR employs a modified Filinov filtration (MFF) scheme to control the overall phase of the SC integrand, extending the applicability of SC methods to complex systems while retaining their ability to accurately describe quantum coherence effects. Here, we address questions regarding the effectiveness of the MFF scheme in combination with SC dynamics. Previous work showed that this filtering scheme is of limited utility in the context of semiclassical wavepacket propagation, but we find the MFF is extraordinarily powerful in the context of correlation functions. By examining trajectory phase and amplitude contributions to the real-time SC correlation function in a model system, we clearly demonstrate that the MFF serves to reduce noise by damping amplitude only in regions of highly oscillatory phase leading to a reduction in computational effort while retaining accuracy. Further, we introduce a novel and efficient MQC-IVR formulation that allows for linear scaling in computational cost with the total simulation length, a significant improvement over the more-than quadratic scaling exhibited by the original method.

Tables

Table 1. Table 1: The total number of trajectories, N t ​ r ​ a ​ j subscript 𝑁 𝑡 𝑟 𝑎 𝑗 N_{traj} , required to achieve numerical convergence for the position correlation function, as computed with different semiclassical formulations.

IVR formulation			$c$			$N_{traj}$
DHK			0			$3.0\times {10}^6$
			$0.7$			$2.4\times {10}^4$
DF-MQC			$3.0$			$9.6\times {10}^3$
			$500.0$			$6.0\times {10}^2$
Husimi			$\infty$			$2.4\times {10}^2$

Table 2. Table 2: A comparison of CPU time to calculate the MQC-IVR position correlation function with the two implementations and c = 0.7 𝑐 0.7 c=0.7 .

Implementation				CPU Time /seconds
DF				$2$
FB				$288$

Equations

$$I=\int \text{d}\mathbf{z}g(\mathbf{z})e^{i\varphi \left(\mathbf{z}\right)}$$

$$I(\mathbf{c})=\int \text{d}\mathbf{z}F(\mathbf{z};\mathbf{c})g(\mathbf{z})e^{i\varphi (\mathbf{z})},$$

$$F(\mathbf{z},\mathbf{c})=\det {\left|\mathbb{I}+i\mathbf{c}\frac{{\partial }^2\varphi }{\partial {\mathbf{z}}^2}\right|}^{\frac{1}{2}}{e}^{-\frac{1}{2}{\frac{\partial \varphi }{\partial \mathbf{z}}}^T\mathbf{c}\frac{\partial \varphi }{\partial \mathbf{z}}},$$

$$\underset{\mathbf{c}\to 0}{\lim }I(\mathbf{c})=I,$$

$$\underset{\mathbf{c}\to \infty }{\lim }I(\mathbf{c})=\sum _{j}g({\mathbf{z}}_j)e^{i\varphi ({\mathbf{z}}_j)}{\left[\det {\left|\frac{1}{2\pi i}\frac{{\partial }^2\varphi }{\partial {\mathbf{z}}^2}\right|}_{\mathbf{z}={\mathbf{z}}_j}\right]}^{-\frac{1}{2}},$$

$$C_{AB}(t)=\text{Tr}\left[\hat{\text{A}}{e}^{\frac{i}{\hslash }\hat{H}t}\hat{\text{B}}{e}^{-\frac{i}{\hslash }\hat{H}t}\right].$$

$$C_{AB}(t)=$$

$$\times ⟨{\mathbf{p}}_0{\mathbf{q}}_0⟩\hat{\text{A}}{\mathbf{p}}_0^{\prime }{\mathbf{q}}_0^{\prime }{e}^{i[S_t({\mathbf{p}}_0,{\mathbf{q}}_0)+S_{-t}({\mathbf{p}}_t^{\prime },{\mathbf{q}}_t^{\prime })]/\hslash }$$

$$\times ⟨{\mathbf{p}}_t^{\prime }{\mathbf{q}}_t^{\prime }⟩\hat{\text{B}}{\mathbf{p}}_t{\mathbf{q}}_t{C}_t({\mathbf{p}}_0,{\mathbf{q}}_0)C_{-t}({\mathbf{p}}_t^{\prime },{\mathbf{q}}_t^{\prime }),$$

$$⟨\mathbf{x}⟩\mathbf{pq}={\left(\frac{\det \left|\gamma \right|}{{\pi }^N}\right)}^{\frac{1}{4}}{e}^{-\frac{1}{2}(\mathbf{x}-\mathbf{q}{)}^T\gamma (\mathbf{x}-\mathbf{q})+\frac{i}{\hslash }{\mathbf{p}}^T(\mathbf{x}-\mathbf{q})},$$

$$C_t^2({\mathbf{p}}_0,{\mathbf{q}}_0)=$$

$$\displaystyle\left.-i\hbar\mathbf{\gamma}_{t}^{\frac{1}{2}}\mathbf{M}_{qp}^{f}\mathbf{\gamma}_{0}^{\frac{1}{2}}+\frac{i}{\hbar}\mathbf{\gamma}_{t}^{-\frac{1}{2}}\mathbf{M}_{pq}^{f}\mathbf{\gamma}_{0}^{-\frac{1}{2}}\right]\bigg{|},$$

$${\mathbf{\Delta }}_{p_t}={\mathbf{p}}_t^{\prime }-{\mathbf{p}}_t\text{and}{\mathbf{\Delta }}_{q_t}={\mathbf{q}}_t^{\prime }-{\mathbf{q}}_t,$$

$$C_{AB}(t)=$$

$$\times ⟨{\mathbf{p}}_0{\mathbf{q}}_0⟩\hat{\text{A}}{\mathbf{p}}_0^{\prime }{\mathbf{q}}_0^{\prime }{e}^{i[S_t({\mathbf{p}}_0,{\mathbf{q}}_0)+S_{-t}({\mathbf{p}}_t^{\prime },{\mathbf{q}}_t^{\prime })]}$$

$$\times ⟨{\mathbf{p}}_t^{\prime }{\mathbf{q}}_t^{\prime }⟩\hat{\text{B}}{\mathbf{p}}_t{\mathbf{q}}_{t}D({\mathbf{p}}_0,{\mathbf{q}}_0,{\mathbf{\Delta }}_{p_t},{\mathbf{\Delta }}_{q_t},{\mathbf{c}}_p,{\mathbf{c}}_q)$$

$$\times e^{-\frac{1}{2}{\mathbf{\Delta }}_{q_t}^T{\mathbf{c}}_q{\mathbf{\Delta }}_{q_t}}{e}^{-\frac{1}{2}{\mathbf{\Delta }}_{p_t}^T{\mathbf{c}}_p{\mathbf{\Delta }}_{p_t}},$$

$$C_{AB}(t)=$$

$$⟨{\mathbf{p}}_0{\mathbf{q}}_0⟩\hat{\text{A}}{\mathbf{p}}_0{\mathbf{q}}_0⟨{\mathbf{p}}_t{\mathbf{q}}_t⟩\hat{\text{B}}{\mathbf{p}}_t{\mathbf{q}}_t,$$

$$V(x)=x^2-0.1x^3+0.1x^4.$$

$$\Psi (x;t=0)={\left(\frac{{\gamma }_x}{\pi }\right)}^{\frac{1}{4}}{e}^{-\frac{{\gamma }_x}{2}(x-q_i{)}^2+i{p}_i(x-q_i)},$$

$$⟨x{⟩}_t=$$

$$\times e^{i[S_t(p_0,q_0)+S_{-t}(p_t^{{}^{\prime }},q_t^{{}^{\prime }})]}⟨p_0{q}_0⟩p_i{q}_i⟨p_i{q}_i⟩p_0^{{}^{\prime }}{q}_0^{{}^{\prime }}$$

$$\times q_t{e}^{-\frac{1}{2}{c}_p{\Delta }_{p_t}^2},$$

$$\omega (p_0,q_0,{\Delta }_{p_t};c_p)=F|⟨p_0{q}_0⟩p_i{q}_i{|}^2{e}^{-\frac{1}{2}{c}_p{\Delta }_{p_t}^2},$$

$$\left|1-E(t)/E(0)\right|$$

$$C_{AB}(t)$$

$$\omega ({\mathbf{p}}_0,{\mathbf{q}}_0,$$

$$\braket{\mathbf{p}^{\prime}_{t}\mathbf{q}^{\prime}_{t}}{\hat{x}}{\mathbf{p}_{t}\mathbf{q}_{t}}=\frac{1}{2}\big{[}(q_{tx}^{\prime}+q_{tx})-\frac{i}{\gamma_{xx}}\Delta_{p_{tx}}\big{]}\braket{\mathbf{p}^{\prime}_{t}\mathbf{q}^{\prime}_{t}}{\mathbf{p}_{t}\mathbf{q}_{t}}.$$

$$D({\mathbf{p}}_0,{\mathbf{q}}_0,{\mathbf{\Delta }}_{p_t},{\mathbf{\Delta }}_{q_t};{\mathbf{c}}_p,{\mathbf{c}}_q)=2^{-\frac{N}{2}}\det ({\gamma }_0^{-1}\mathbf{G}{)}^{\frac{1}{2}}$$

$$\displaystyle\times\det\bigg{[}\frac{1}{2}(\mathbf{M}_{pp}^{b}-i\mathbf{\gamma}_{0}\mathbf{M}_{qp}^{b})(\mathbf{G}^{-1}+\mathbb{I})(\mathbf{M}_{pp}^{f}\mathbf{\gamma}_{0}+i\mathbf{M}_{pq}^{f})$$

$$+({\gamma }_0{\mathbf{M}}_{qq}^b+i{\mathbf{M}}_{pq}^b)(\frac{1}{2}{\gamma }_t^{-1}+{\mathbf{c}}_p){\mathbf{G}}^{-1}({\mathbf{M}}_{pp}^f{\gamma }_0+i{\mathbf{M}}_{pq}^f)$$

$$+\frac{1}{2}({\gamma }_0{\mathbf{M}}_{qq}^b+i{\mathbf{M}}_{pq}^b)({\mathbf{G}}^{-1}+\mathbb{I})({\mathbf{M}}_{qq}^f-i{\mathbf{M}}_{qp}^f{\gamma }_0)$$

$$\displaystyle+(\mathbf{M}_{pp}^{b}-i\mathbf{\gamma}_{0}\mathbf{M}_{qp}^{b})(\frac{1}{2}\mathbf{\gamma}_{t}+\mathbf{c}_{q})\mathbf{G}^{-1}(\mathbf{M}_{qq}^{f}-i\mathbf{M}_{qp}^{f}\mathbf{\gamma}_{0})\bigg{]}^{\frac{1}{2}},$$

$$D_q({\mathbf{p}}_0,{\mathbf{q}}_0,{\mathbf{\Delta }}_{p_t},{\mathbf{c}}_p)=$$

$$\displaystyle\det\bigg{|}\mathbf{\gamma}_{0}^{-1}\mathbf{c}_{p}\bigg{[}(\mathbf{M}_{pp}^{b}-i\mathbf{\gamma}_{0}\mathbf{M}_{qp}^{b})(\mathbf{M}_{pp}^{f}\mathbf{\gamma}_{0}+i\mathbf{M}_{pq}^{f})$$

$$+({\gamma }_0{\mathbf{M}}_{qq}^b+i{\mathbf{M}}_{pq}^b)({\mathbf{M}}_{qq}^f-i{\mathbf{M}}_{qp}^f{\gamma }_0)$$

$$\displaystyle+\;(\mathbf{M}_{pp}^{b}-i\mathbf{\gamma}_{0}\mathbf{M}_{qp}^{b})\mathbf{c}_{p}^{-1}(\mathbf{M}_{qq}^{f}-i\mathbf{M}_{qp}^{f}\mathbf{\gamma}_{0})\bigg{]}\bigg{|}^{\frac{1}{2}}$$

Full text

††preprint: AIP/123-QED

I Introduction

Semiclassical (SC) approximations to exact, real-time quantum correlation functions such as the Double Herman-Kluk Initial Value Representation hk1 (DHK-IVR) are unique in their ability to accurately describe a wide range of chemical processes rev1 ; rev2 ; rev3 ; rev4 including systems where nuclear and electronic quantum coherences play a significant dynamic role coh ; dhk2 such as in the vibrational motion following a Franck-Condon excitation, franck or in various nonadiabatic processes. ananth ; nonad1 ; nonad2 ; nonad3 ; nonad4 Unfortunately, the computational cost of converging the highly oscillatory SC integrand scales exponentially with system dimensionality, limiting the application of these methods to relatively low-dimensional systems.

Efforts to mitigate the so-called SC ‘sign’ problem have led to the development of several forward-backward techniques fb1 ; fb2 ; fb3 ; fb4 ; fb5 ; fb6 ; fb7 that offer a small reduction in computational cost while retaining accuracy. However, the simulation of truly complex systems remains limited to methods that do not fully describe coherence. This includes SC methods that employ an additional linearization lsc1 ; lsc2 or related approximations related1 ; related2 ; fbsd1 ; fbsd2 , as well as model path-integral based dynamics such as Centroid Molecular Dynamics, cmd1 ; cmd2 and Ring Polymer Molecular Dynamics and its extensions. rpmd1 ; rpmd2 ; rpmd3 ; rpmd4 ; rpmd5 ; rpmd6 More recently, a handful of exact qcpi and approximate methods that can account for electronic coherence effects mvrpmd ; jess1 ; cotton ; tao_cc at short times have emerged, however the inclusion of nuclear coherences in high-dimensional systems remains out of reach.

To address this challenge, two of us recently introduced a new SC method for the calculation of real-time correlation functions: the Mixed Quantum-Classical Initial Value Representation (MQC-IVR). mqcivr This method is derived using the modified Filinov filtration (MFF) technique filinov1 ; filinov2 ; filinov3 ; filinov4 ; filinov5 ; filinov6 ; filinov7 ; filinov8 ; filinov9 ; filinov10 to associate a filtration parameter with each mode of the system; modulating these parameters allows us to control the strength of the filter and hence the extent of quantization of an individual mode. As demonstrated in our introductory paper, for an optimal choice of parameters, MQC-IVR correlation functions correctly account for quantum coherence effects and require significantly fewer trajectories to achieve numerical convergence than standard quantum-limit SC methods like the DHK-IVR.

In this paper, we analyze the effectiveness of the MFF technique, addressing concerns raised by an earlier study that found the MFF to be of limited utility in the context of SC wavepacket propagation. batista Specifically, the study found the MFF to be ineffective for both low-dimensional non-chaotic systems, where the amplitude of the SC integrand in regions of highly oscillatory phase was shown to be negligible, and for chaotic systems where the MFF was necessary to achieve numerical convergence but led to loss of valuable phase information. batista Here, we undertake a similar analysis of the MQC-IVR correlation function where phase contributions arise from the net action of a pair of forward-backward trajectories rather than the action along a single trajectory. We find that in a simple low-dimensional system, (i) the trajectory pairs with highly oscillatory phase make numerical convergence difficult while contributing little to the actual ensemble average, and (ii) the MFF technique successfully filters contributions from pairs of trajectories that contribute highly oscillatory terms to the MQC-IVR integrand, significantly reducing computational cost while retaining accuracy. Having validated our filtration scheme, we further reformulate the MQC-IVR expression to facilitate a novel, highly parallel double-forward (DF) implementation that exhibits linear scaling in CPU time with simulation length, unlike the more-than quadratic scaling exhibited by our original forward-backward (FB) implementation.

The paper is organized as follows: in Section II, we begin with a brief review of the MFF technique and the original MQC-IVR formalism. In Section III, we introduce a simple 1D model, provide implementation details, and present our analysis of the phase and amplitude of the MQC-IVR integrand, demonstrating the effectiveness of the MFF scheme. In Section IV, we introduce our modified MQC-IVR expression with the novel DF implementation and numerically calculate the position correlation function, showing a two-order of magnitude reduction in CPU time for our 1D model system. We conclude in Section V.

II Theory

II.1 Modified Filinov Filtration

The MFF technique filinov1 ; filinov2 ; filinov3 can be used to smooth oscillatory integrals of the form,

[TABLE]

where the functions $g(\mathbf{z})$ and $\phi(\mathbf{z})$ are complex-valued. Using the MFF to smooth the oscillatory integrand, we obtain

[TABLE]

where the prefactor is defined as

[TABLE]

and $\mathbb{I}$ is the identity matrix. In Eq. (2), the filtering strength is determined by elements of the diagonal matrix of tuning parameters, $c_{ii}\geq 0$. In the limit of very small tuning parameters the integrand is unchanged,

[TABLE]

and in the limit of very large tuning parameter values the MFF is equivalent to a stationary phase approximation,

[TABLE]

where the subscript $j$ indexes stationary phase points. Modifying the values of $c_{ii}$ thus controls the extent to which the $i^{\text{th}}$ system mode is filtered, and can, ideally, be chosen to facilitate fast numerical convergence of the integral without loss of valuable phase information.

II.2 Forward-backward MQC-IVR

The quantum real-time correlation function between two operators $\hat{\text{A}}$ and $\hat{\text{B}}$ for a system with Hamiltonian $\hat{H}$ is written as

[TABLE]

The corresponding semiclassical DHK-IVR approximation to Eq. (6) is a double phase space average obtained by replacing both the forward and backward time propagators with HK-IVR propagators, hk1 ; filinov10 ; rev3 ; rev4

[TABLE]

where $N$ is the system dimensionality, the forward and backward trajectories have initial phase space points, $(\mathbf{p}_{0},\mathbf{q}_{0})$ and $(\mathbf{p}^{\prime}_{t},\mathbf{q}^{\prime}_{t})$, final phase space points, $(\mathbf{p}_{t},\mathbf{q}_{t})$ and $(\mathbf{p}^{\prime}_{0},\mathbf{q}^{\prime}_{0})$, and classical actions, $S_{t}$ and $S_{-t}$, respectively. The position representation of the coherent states employed in Eq. (II.2) is given by

[TABLE]

where the diagonal $N\times N$ matrix, $\mathbf{\gamma}$, determines the spread of the coherent state in phase space. Finally, in Eq. (II.2), the prefactor for the forward trajectory is

[TABLE]

where $\mathbf{M}_{\alpha\beta}^{f}=\frac{\partial\mathbf{\alpha}_{t}}{\partial\mathbf{\beta}_{0}}$ are the forward monodromy matrices, and the prefactor for the backward trajectory, $C_{-t}$, is similarly defined in terms of the backward monodromy matrices, $\mathbf{M}_{\alpha\beta}^{b}=\frac{\partial\mathbf{\alpha}_{0}^{\prime}}{\partial\mathbf{\beta}_{t}^{\prime}}$.

Non-zero phase space displacements between a pair of forward and backward trajectories at time $t$,

[TABLE]

leads to a net non-zero action, $S_{t}+S_{-t}$, that contributes to the overall phase of a trajectory pair. Pairs of trajectories with zero phase space displacement have zero net action and make no contribution to the overall phase while widely diverging trajectory pairs make the semiclassical integrand in Eq. (II.2) highly oscillatory and difficult to numerically converge. This motivates the use of a filtering technique like the MFF to smooth the oscillatory integrand by controlling the phase space displacement between forward and backward trajectories.

Applying the MFF to the oscillatory DHK-IVR integrand in Eq. (II.2) yields the MQC-IVR approximation for a real-time correlation function mqcivr

[TABLE]

where we set $\hbar=1$ here and in the rest of the manuscript. In Eq. (II.2), the prefactor $D(\mathbf{p}_{0},\mathbf{q}_{0},\mathbf{\Delta}_{p_{t}},\mathbf{\Delta}_{q_{t}},\mathbf{c}_{p},\mathbf{c}_{q})$ is a determinant of an $N\times N$ matrix containing elements of the forward and backward monodromy matrices (full definition provided in Appendix A.1), and $\mathbf{c}_{p}$ and $\mathbf{c}_{q}$ are $N\times N$ diagonal matrices of tuning parameters. Matrix elements of the tuning parameter $\mathbf{c}_{p}$ control the extent of momentum displacement in a particular mode between trajectory pairs, and the elements of $\mathbf{c}_{q}$ similarly control position displacement. In the limit $\mathbf{c}_{p},\mathbf{c}_{q}\rightarrow 0$, the MQC-IVR expression in Eq. (II.2) becomes identical to the DHK-IVR expression in Eq. (II.2). As $\mathbf{c}_{p},\mathbf{c}_{q}\rightarrow\infty$, the path displacements in Eq. (10) become zero as does the overall action of the forward-backward trajectory pairs. In this limit, the MQC-IVR correlation function reduces to the Husimi-IVR,

[TABLE]

a classical-limit IVR method that employs Husimi functions hus1 of the two operators rather than the Wigner functions used in the LSC-IVR formulation. The Husimi-IVR correlation function involves a classical average, making it an appealing method for high-dimensional systems but fails to account for quantum effects like deep tunneling and coherence. Choosing intermediate values for the matrix elements of $\mathbf{c}_{p}$ and $\mathbf{c}_{q}$ thus allows us to control the extent of ‘quantum’-ness incorporated in our simulations, and choosing different values for different modes allows us to perform mode-specific quantization within a uniform semiclassical dynamic framework.

III Analyzing the MFF

III.1 Model System

For clarity, we analyze the effectiveness of the MFF using the 1D anharmonic oscillator employed as a model system in our introductory MQC-IVR paper, mqcivr

[TABLE]

The initial state is chosen to be a coherent state, $\ket{p_{i}q_{i}}$, and the corresponding wavefunction is given by

[TABLE]

with $q_{i}=1$, $p_{i}=0$, and $\gamma_{x}=\sqrt{2}$, all in atomic units.

III.2 FB Implementation

We calculate the time-dependent position expectation value with $\hat{\text{A}}=\ket{p_{i}q_{i}}\bra{p_{i}q_{i}}$, the projection operator corresponding to an initial coherent state, and $\hat{\text{B}}=\hat{x}$, the position operator. In the MQC-IVR formulation, this expectation value is written as mqcivr

[TABLE]

where the prefactor $D_{q}(p_{0},q_{0},\Delta_{p_{t}},c_{p})$ is provided in Appendix A.2. Note that the integral over $\Delta_{q_{t}}$ has been evaluated analytically in Eq. (III.2) using the $\gamma_{t}\rightarrow\infty$ limit. mqcivr

Initial positions, momenta, and the momentum jump are sampled from the distribution,

[TABLE]

with $F=\sqrt{\frac{c_{p}}{(2\pi)^{3}}}$ for normalization, and time-evolved using a fourth-order symplectic integrator integrator with a time step of $\Delta t\,=\,0.05\,\text{a.u.}$ along with the monodromy matrices. Energy conservation along all trajectories is monitored with a tolerance parameter, $\epsilon=10^{-4}$, such that

[TABLE]

We also track the Maslov index of the prefactor in order to select the correct branch of the complex square root.

III.3 Phase and Amplitude Analysis

Since the position correlation function defined in Eq. (III.2) involves only a momentum displacement, $\Delta_{p_{t}}$, we can analyze the effectiveness of the MFF by examining the phase and amplitude contributions to the MQC-IVR integrand from different forward-backward trajectory pairs that share the same initial conditions at time $t=0$, and differ only by the $\Delta_{p_{t}}$ values. We choose initial conditions $q_{0}=1$ a.u. and $p_{0}=0$ a.u. corresponding to the center of the initial coherent state, propagate forward for time $t$, introduce a change in momentum, $\Delta_{p_{t}}$, and then propagate backward to time $t=0$. We plot the amplitude and phase of the integrand for each pair of forward-backward trajectories as a function of $\Delta_{p_{t}}$ in Fig. 1(a) and Fig. 1(b) for $t=22$ a.u. and in Fig. 1(c) and Fig. 1(d) for $t=61\,\text{a.u.}$. The two times chosen here correspond to values at which the classical and quantum position correlation functions have significantly different amplitudes.

The different colored lines in the phase and amplitude plots in Fig. 1 correspond to different values of the Filinov parameter, $c_{p}$, ranging from the quantum-limit (shown in purple) to the classical-limit (shown in red). Consider the phase corresponding to the momentum jump value $\Delta_{p_{t}}=0$ in Fig. 1(a) and Fig. 1(c). At this point, the forward and backward trajectories coincide and the net phase is zero (or $\pm\pi$ due to a phase contribution from operator $\hat{\text{B}}$) – these are the trajectory pairs that contribute to the classical-limit Husimi-IVR correlation function in Eq. (II.2), as indicated by the narrow red ($c_{p}=200$) peak in the corresponding amplitude plots in Fig. 1(b) and Fig. 1(d). As the magnitude of $\Delta_{p_{t}}$ increases, the path of the forward and backward trajectories are significantly different, resulting in a net non-zero phase contribution that varies rapidly with small changes in momentum jump. Although the regions of slowly varying phase in Fig. 1(a) and Fig. 1(c) (enclosed within solid black lines) are relatively invariant to the choice of $c_{p}$ values, as we go from low $c_{p}$ values that describe quantum-limit behavior to high $c_{p}$ values for classical-limit behavior, the amplitude changes significantly, approaching zero for large values of $\Delta_{p_{t}}$. Our analysis of the MFF based on a single trajectory pair thus suggests that we can eliminate contributions from regions of highly oscillatory phase, $|\Delta_{p_{t}}|\geq 3$ in both Fig. 1(a) and Fig. 1(c).

Next we analyze the amplitude and phase of an average integrand obtained by summing over contributions from an ensemble of forward-backward trajectory pairs to verify that as we approach numerical convergence, regions of highly oscillatory phase do, in fact, have a negligible contribution to the ensemble average correlation function. Using a weak filter (quantum-limit, $c_{p}=0.05$) and simulation time $t=22$ a.u., we plot the average phase of the MQC-IVR integrand in Fig. 2(a) as a function of the momentum jump. In Fig. 2(a) we find a region of slowly varying phase corresponding to the region $-4$ a.u.$\leq\Delta_{p_{t}}\leq 4$ a.u (enclosed between two vertical black lines). The integrand amplitude is plotted in Fig 2(b); we find that as the number of trajectories included in the ensemble average is increased, the amplitude vanishes in regions of highly oscillatory phase. Taken together, these observations suggest that trajectory pairs with $\Delta_{p_{t}}$ values outside the slowly varying phase zone ($|\Delta_{p_{t}}|\geq 4$ a.u.) contribute very little to the overall correlation function and make it difficult to achieve numerical convergence. This makes a strong case for the use of MFF to increase computational efficiency – for instance, with a stronger filter, say $c_{p}=0.7$ (green Gaussian in Fig. 2(b)), it becomes possible to explicitly exclude trajectory pairs that contribute only noise by making the amplitude zero. We show, in Fig. 3(b) of the following section, that an MQC-IVR simulation employing $c_{p}=0.7$ does indeed recover quantum recurrence in the position correlation function at long times with significantly fewer trajectories.

IV The Double-Forward MQC-IVR

IV.1 Re-formulating MQC-IVR

In the original forward-backward (FB) implementation, mqcivr with general operators $\hat{\text{A}}$ and $\hat{\text{B}}$, we sample the initial phase space points of the forward trajectory $(\mathbf{p}_{0},\mathbf{q}_{0})$ and the path displacement variables at time $t$, ($\mathbf{\Delta}_{p_{t}},\mathbf{\Delta}_{q_{t}}$), to evaluate the integral in Eq. (II.2). To calculate the MQC-IVR correlation at $N_{t}$ time points requires propagating forward trajectories for $N_{t}$ time steps, using the sampled phase space jump values to generate new initial conditions for the backward trajectories at each time step, and then propagating trajectories backward for $N_{t}$ steps. Thus, a correlation function with $N_{t}$ time points involves a total of $\frac{1}{2}(N_{t}^{2}+3N_{t})$ propagation steps for each trajectory pair.

Here, we introduce a novel double-forward (DF) implementation that significantly reduces the computational effort involved in an MQC-IVR calculation. We re-write the MQC-IVR correlation function in Eq. (II.2) such that the variables of integration are the initial phase space points of two independent forward trajectories,

[TABLE]

Details of the derivation are provided in Appendix B.1 along with the explicit form of the prefactor, $D_{t}$. As expected, the limits of the MQC-IVR correlation function remain the same as before yielding the DHK-IVR when $\mathbf{c}_{p},\mathbf{c}_{q}\rightarrow 0$ and the Husimi-IVR when $\mathbf{c}_{p},\mathbf{c}_{q}\rightarrow\infty$. We evaluate the correlation function defined in Eq. (18) by sampling the variables $(\mathbf{p}_{0},\mathbf{q}_{0})$ and $(\mathbf{p}^{\prime}_{0},\mathbf{q}^{\prime}_{0})$ from a correlated sampling distribution,tao

[TABLE]

where we define mean variables, $\mathbf{\bar{x}}=\frac{\mathbf{x}_{0}^{\prime}+\mathbf{x}_{0}}{2}$, and we introduce trajectory displacement variables at time $t=0$, $\mathbf{\Delta}_{p_{0}}=\mathbf{p}^{\prime}_{0}-\mathbf{p}_{0}$ and $\mathbf{\Delta}_{q_{0}}=\mathbf{q}^{\prime}_{0}-\mathbf{q}_{0}$. For normalization we have $F=\frac{\sqrt{\det|\mathbf{c}_{p}|\det|\mathbf{c}_{q}|}}{(4\pi^{2})^{N}}$. Calculating the MQC-IVR correlation function now involves propagating two independent trajectories forward in time. The contribution to a correlation function of length $N_{t}$ is obtained simply from a given pair of forward trajectories using a total of $2N_{t}$ propagation steps.

Thus, the DF implementation exhibits linear scaling in CPU time with total simulation length relative to the more-than quadratic scaling of the original FB formulation. Furthermore, the use of two independently propagated forward trajectories, rather than a forward-backward structure, allows for highly parallel implementation.

IV.2 Simulation Details

Initial conditions for the forward trajectory pairs are sampled using the correlated sampling function in Eq. (IV.1) with $N=1$. As before, we calculate the position correlation function, $\langle x\rangle_{t}$, where $\hat{\text{A}}=\ket{p_{i}q_{i}}\bra{p_{i}q_{i}}$ and $\hat{\text{B}}=\hat{x}$ using the DF implementation in Eq. (18). The matrix element of the position operator is given by

[TABLE]

Details of trajectory integration remain the same as before. The SC-Corr code package, corrcode developed in-house and available as open-source software, was used to perform these calculations.

IV.3 Numerical Results for DF Correlation Functions

In Fig. 3(a), we plot $\langle x\rangle_{t}$ as a function of time obtained from the DF implementation with different values of the tuning parameters $c_{q}=c_{p}=c$. Comparing this result against exact quantum, DHK-IVR, and Husimi-IVR simulations, we demonstrate that the new DF implementation, like the original FB implementation, mqcivr recovers the amplitude recurrences at long times that distinguish the quantum result from the classical result. As the filter strength ($c$ value) is increased, we obtain classical-limit results where the amplitude is damped at long times. In Table 1, we report the number of trajectories required to converge each simulation such that the maximum absolute statistical error across all time steps is $\leq 5\%$. Even for a quantum-limit filter strength, $c=0.7$, the DF implementation captures quantum behavior with two orders of magnitude fewer trajectories than required by a full semiclassical DHK-IVR calculation. We note that the number of trajectories required to achieve numerical convergence with the DF implementation is comparable to the original FB implementation, but the CPU time required per trajectory pair is significantly lower.

An additional advantage of the DF implementation becomes apparent in Fig. 3(b) and Fig. 3(c) where we compare the position correlation function obtained from the DF and FB MQC-IVR implementations for $c=0.7$ and $c=10$ respectively. The DF formulation better captures quantum recurrences at long times since the filter now determines phase space displacements at time zero, allowing a broader range of trajectory displacements to be included at time $t$. Further, the total CPU time required to compute $\langle x\rangle_{t}$ with $c=0.7$, $2.4\times 10^{4}$ trajectories, and $1600$ time steps for both implementations is reported in Table 2, demonstrating that the DF implementation increases computational efficiency by two orders of magnitude.

V Conclusions

We show that the MFF technique used in MQC-IVR acts to apriori filter out phase contributions from widely diverging trajectory pairs that contribute little, on average, to the correlation function. This makes it an effective and efficient tool for enhancing numerical convergence in the calculation of correlation functions without loss of important quantum information.

Further, we introduce a novel DF implementation for the MQC-IVR correlation function that scales linearly with simulation length compared to the quadratic scaling of the FB implementation. This, combined with an improved ability to capture quantum recurrences at long times, immediately extends the applicability of MQC-IVR to high-dimensional system simulations.

Acknowledgements

This work is supported in part by the Army Research Office Grant No. W911NFD-13-1-0102 and in part by Cornell University Start-Up funding. N. A. additionally acknowledges support from the Research Corporation for Science through a Cottrell Scholar Award and a Sloan Foundation Research Fellowship.

Appendix A Forward-backward MQC-IVR prefactors

A.1 General operator $\hat{\text{B}}$

The MQC-IVR prefactor for a general operator $\hat{\text{B}}$ is given by,

[TABLE]

with diagonal matrix $\mathbf{G}=(\mathbf{c}_{q}+\mathbf{\gamma}_{t})\mathbf{c}_{p}+\mathbf{c}_{q}(\mathbf{\gamma}_{t}^{-1}+\mathbf{c}_{p})$. The superscript $f$ denotes the monodromy matrix element corresponding to the forward trajectory beginning at $(\mathbf{p}_{0},\mathbf{q}_{0})$, while the superscript $b$ denotes the monodromy matrix element corresponding to the backward trajectory beginning at $(\mathbf{p}^{\prime}_{t},\mathbf{q}^{\prime}_{t})$.

A.2 Position-space operator $\hat{\text{B}}=\hat{\text{B}}(\hat{q})$

The MQC-IVR prefactor for a position-space operator $\hat{\text{B}}(\hat{q})$ is obtained by evaluating the $\gamma_{t}\rightarrow\infty$ limit,

[TABLE]

Appendix B Double Forward MQC-IVR Formulation

B.1 Changing variables of integration

We start with the DHK-IVR of Eq. (II.2) and change variables of integration from $(\mathbf{p}^{\prime}_{t},\mathbf{q}^{\prime}_{t})$ to $(\mathbf{p}^{\prime}_{0},\mathbf{q}^{\prime}_{0})$ using Liouville’s theorem. The action of the backward trajectory, $S_{-t}(\mathbf{p}^{\prime}_{t},\mathbf{q}^{\prime}_{t})$, is then replaced by its forward counterpart, $S_{t}(\mathbf{p}^{\prime}_{0},\mathbf{q}^{\prime}_{0})$, and we use the following monodromy matrix identity for the backward trajectory in order to obtain a new prefactor expression,

[TABLE]

with $\mathbf{M}_{\alpha\beta}^{f\prime}=\frac{\partial\alpha_{t}^{\prime}}{\partial\beta_{0}^{\prime}}$ The resulting expression for the correlation function is a double-forward DHK-IVR,

[TABLE]

We now follow the original MQC-IVR derivation, mqcivr implementing the MFF scheme to smooth the oscillatory integrand in Eq. (21). We take $\phi$ to include the coherent state exponentials in Eq. (21) and the action terms, excluding any contributions from operators $\hat{\text{A}}$ and $\hat{\text{B}}$ as well as any phase contributions from the prefactors,

[TABLE]

The derivation proceeds much like the original method (from Eq. (12) onward in the original manuscript mqcivr ), but the vector of phase space displacements is now given by

[TABLE]

and the block diagonal matrix of tuning parameters is

[TABLE]

After simplifying the integrand we obtain Eq. (18) and the prefactor takes the form,

[TABLE]

with diagonal matrix $\mathbf{G}=(\mathbf{c}_{q}+\mathbf{\gamma}_{0})\mathbf{c}_{p}+\mathbf{c}_{q}(\mathbf{\gamma}_{0}^{-1}+\mathbf{c}_{p})$. We note that unlike the FB implementation, the DF implementation will typically always include both position and momentum jumps even if operator $\hat{\text{B}}$ is a pure position or momentum operator.