Semiclassical dynamics in the mixed quantum-classical limit

TL;DR

This paper introduces the AMQC-IVR method, a mixed quantum-classical semiclassical approach that improves computational efficiency and accuracy in calculating quantum correlation functions and reaction rates, especially in weakly coupled systems.

Contribution

It analytically derives the AMQC-IVR expression, combining filtered classical modes with full semiclassical quantum modes, and demonstrates its effectiveness through numerical examples.

Findings

AMQC-IVR accurately computes quantum correlation functions.

The method is efficient for systems with weak quantum-classical coupling.

A separable prefactor approximation reduces computational cost in specific regimes.

Abstract

The semiclassical Double Herman-Kluk Initial Value Representation is an accurate approach to computing quantum real time correlation functions, but its applications are limited by the need to evaluate an oscillatory integral. In previous work, we have shown that this `sign problem' can be mitigated using the modified Filinov filtration technique to control the extent to which individual modes of the system contribute to the overall phase of the integrand. Here we follow this idea to a logical conclusion: we analytically derive a general expression for the mixed quantum-classical limit of the semiclassical correlation function - AMQC-IVR, where the phase contributions from the `classical' modes of the system are filtered while the `quantum' modes are treated in the full semiclassical limit. We numerically demonstrate the accuracy and efficiency of the AMQC-IVR formulation in calculations of quantum correlation functions and reaction rates using three model systems with varied coupling strengths between the classical and quantum subsystems. We also introduce a separable prefactor approximation that further reduces the computational cost, but is only accurate in the limit of weak coupling between the quantum and classical subsystems.

Tables

Table 1. Table 1: The time-averaged (relative) % error of each AMQC-IVR result with model 1. The result in parentheses was obtained after quantizing both the anharmonic and harmonic dofs.

$k$	$0.5$	$1.5$	$2.0$
Avg. % Error	$1.07$	$0.96$	$0.91\left(0.31\right)$

Table 2. Table 2: The time-averaged (relative) % error of the SP approximation with model 1.

$k$	$0.5$	$1.5$	$2.0$
Avg. % Error	$1.06$	$1.53$	$1.85$

Table 3. Table 3: The percent error between the AMQC-IVR and SP results for the thermal transmission coefficient at three different coupling strengths.

$\eta /m{\omega }_b$	$0.05$	$0.1$	$0.5$
$\%\text{Error}$	$3.0$	$38$	$65$

Table 4. Table 4: The AMQC-IVR results for the thermal transmission coefficient of model 3.

$\eta /m{\omega }_b$	$\kappa$
0.05	0.58(7)
0.1	0.86(2)
0.2	1.21(4)
0.25	1.31(4)
0.35	1.51(6)
0.4	1.56(2)
0.5	1.62(3)
0.6	1.62(5)
0.75	1.60(2)
0.9	1.54(1)
1.0	1.48(2)
1.1	1.41(1)
1.25	1.30(2)
1.4	1.20(1)
1.5	1.12(1)

Equations

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

$$\times$$

$${\mathbf{z}}_t=$$

$$=$$

$${\mathbf{z}}_t^{\prime }=$$

$$=$$

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

$${\Omega }_{\mathbf{z}{\mathbf{z}}^{\prime }}=⟨\mathbf{z}⟩\hat{\Omega }{\mathbf{z}}^{\prime }.$$

$$\mathbf{c}=\left(\begin{array}{cc}{\mathbf{c}}_{\mathbf{p}}& \mathbb{O}\\ \mathbb{O}& {\mathbf{c}}_{\mathbf{q}}\end{array}\right),$$

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

$$\times$$

$$\displaystyle C_{t}\left(\mathbf{z}_{0}\right)=\det\bigg{[}$$

$$\displaystyle-i\bm{\gamma}_{0}^{\frac{1}{2}}\mathbf{M}_{qp}\bm{\gamma}_{t}^{\frac{1}{2}}+i\bm{\gamma}_{0}^{-\frac{1}{2}}\mathbf{M}_{pq}\bm{\gamma}_{t}^{-\frac{1}{2}}\big{)}\bigg{]}^{\frac{1}{2}},$$

$$C_{AB}(t)=\frac{1}{{\left(2\pi \right)}^N}\int d{\mathbf{z}}_0{A}_{{\mathbf{z}}_0{\mathbf{z}}_0}{B}_{{\mathbf{z}}_t{\mathbf{z}}_t},$$

$$G_t\left({\mathbf{z}}_0,{\mathbf{z}}_0^{\prime };\mathbf{c}\right)=$$

$${\mathbf{z}}_Q=$$

$${\mathbf{z}}_Q^{\prime }=$$

$${\mathbf{z}}_C=$$

$${\mathbf{z}}_C^{\prime }=$$

$$\delta (x)=\underset{a\to \infty }{\lim }{\left(\frac{a}{2\pi }\right)}^{\frac{1}{2}}{e}^{-\frac{a}{2}{x}^2},$$

$$\underset{\genfrac{}{}{0ex}{}{c_{\text{quantum}}\to 0}{c_{\text{classical}}\to \infty }}{\lim }{G}_t\left({\mathbf{z}}_0,{\mathbf{z}}_0^{\prime };\mathbf{c}\right)=$$

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

$$\times$$

$$G_t\left({\mathbf{z}}_0,{\mathbf{z}}_0^{\prime };\mathbf{c}\right)$$

$$G_t\left({\mathbf{z}}_Q,{\mathbf{z}}_Q^{\prime };{\mathbf{c}}_Q\right)e^{i\left[S_t\left({\mathbf{z}}_Q\right)-S_t\left({\mathbf{z}}_Q^{\prime }\right)\right]}$$

$$\times$$

$$C_t\left({\mathbf{z}}_0\right)C_t^{\ast }\left({\mathbf{z}}_0^{\prime }\right)$$

$$C_t\left({\mathbf{z}}_Q\right)C_t^{\ast }\left({\mathbf{z}}_Q^{\prime }\right)e^{i\left[S_t\left({\mathbf{z}}_Q\right)-S_t\left({\mathbf{z}}_Q^{\prime }\right)\right]}$$

$$\times$$

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

$$\times$$

$${\Lambda }_t\left({\mathbf{z}}_0,{\mathbf{z}}_Q^{\prime }\right)=1.$$

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

$$\times$$

$$\hat{H}=$$

$$+\frac{1}{2}{m}_2{\omega }_2^2{\hat{q}}_2^2+k{\hat{q}}_1{\hat{q}}_2,$$

$$\Psi \left(\mathbf{x},0\right)=$$

Full text

††preprint: AIP/123-QED

I Introduction

Semiclassical methods based on the initial value representation (SC-IVR) can be used to characterize quantum mechanical effects of many-body systems in real time.Miller (2001); Thoss and Wang (2004); Kay (2005); Herman and Kluk (1984); Kay (1994a, b, c) Using classical trajectories from molecular dynamics simulations, SC-IVR methods accurately describe bound-state motion, tunneling processes, chemical reaction rates, and coherence effects in both adiabatic and nonadiabatic systems.Venkataraman and Miller (2007); Kay (1997); Sun and Miller (1998); Zhang and Pollak (2004a); Buchholz (2018); Zhang and Pollak (2004b); Moix and Pollak (2009); Elran and Kay (2002); Skinner and Miller (1999); Stock and Thoss (1997); Coronado, Batista, and Miller (2000); Ananth, Venkataraman, and Miller (2007); Miller (2012) And while there are other classes of trajectory-based methods such as ring polymer molecular dynamicsCraig and Manolopoulos (2004, 2005); Menzeleev, Ananth, and Miller III (2011); Habershon et al. (2013); Menzeleev, Bell, and Miller III (2014) and centroid molecular dynamicsCao and Voth (1994); Jang and Voth (1999) that can capture some quantum effects in condensed phase systems, these methods cannot be used for systems where quantum coherence effects play a role. Yet other classes of methods are derived from either the exact path integral representation of the propagator, Makri (1995) the quantum Liouvillian, Kapral and Ciccotti (1999) or wavepacket dynamics, Beck et al. (2000); Meyer, Gatti, and Worth (2009) through a series of rigorous approximations. However, these approaches remain limited to low-dimensional systems or condensed phase systems where a large number of near-classical modes serve to mitigate the importance of long-lived quantum coherence effects. Kapral (2006); Makri (2015); Walters and Makri (2015); Wang and Thoss (2003)

Efforts to make SC-IVR methods computationally feasible focus on dealing with the ‘sign problem’ that arises from the inclusion of a phase from individual trajectories that must be averaged over. There are a number of existing approximations that make SC-IVR theory more amenable to large-scale simulation, such as the widely-used linearized SC-IVR,Wang, Sun, and Miller (1998); Sun, Wang, and Miller (1998); Liu (2015); Shi and Geva (2003) a classical limit of SC theory that is accurate on short time scales but suffers from zero-point energy leakage and fails to describe long-time coherence effects.Habershon and Manolopoulos (2009); Buchholz (2018); Miller (2001); Gelavert et al. (2001); Ananth, Venkataraman, and Miller (2007) Other methods include the various forward-backwardSun and Miller (1999); Wang, Thoss, and Miller (2000); Wang et al. (2001); Gelavert et al. (2001); Thoss, Wang, and Miller (2001) SC-IVRs, and there are promising semiclassical methods based on time-averaging,Kaledin and Miller (2003); Buchholz, Grossmann, and Ceotto (2018) the “divide and conquer” methodology,Ceotto, Di Liberto, and Conte (2017); Di Liberto, Conte, and Ceotto (2018a, b) linearization,Lee, Huo, and Coker (2016) symmetrical windowing,Cotton and Miller (2013) and semiclassical quantizationLoring (2017) for the calculation of electronic coherence and/or 1D and 2D vibrational and vibronic spectra. Still, however, there is need of practical methods that include a true description of nuclear coherence in real time for the study of interesting processes such as the generation of hot-electrons at metal surfaces, intramolecular vibrational relaxation, molecular collisions, and other electronically and vibronically nonadiabatic processes.Uzer and Miller (1991); Golibrzuch et al. (2015); Kruger et al. (2015); Domcke and Yarkony (2012); Reese and Nocera (2009); Gray and Winkler (1996); Hammes-Schiffer and Stuchebrukhov (2010)

The mixed quantum-classical IVR (MQC-IVR) is a semiclassical approach to calculating correlation functions that has shown promise in mitigating the SC-IVR sign problem in low-dimensional adiabatic and nonadiabatic systems.Antipov, Ye, and Ananth (2015); Church, Antipov, and Ananth (2017); Church et al. (2018) It is derived by applying modified Filinov filtrationFilinov (1986); Makri and Miller (1987); Brewer, Hulme, and Manolopoulos (1997); Wang, Manolopoulos, and Miller (2001); Spanner, Batista, and Brumer (2005) (MFF) to the double Herman-Kluk (DHK-IVR) formulation of the correlation function, making the level of theory used on each dof tunable (via an adjustable ‘tuning’ parameter) between classical and quantum limits of SC-IVR theory.222Throughout the study we may refer to dofs treated in the classical limit of SC-IVR theory as ‘classical’, and to dofs treated in the quantum limit of SC-IVR theory as ‘quantum’ in order to distinguish the level of SC-IVR theory used for different modes in a given system. We do this with the understanding that the whole system is described at the semiclassical level of theory. The zero-limit of the tuning parameters is equivalent to treating the full system at the DHK-IVR level of theory (i.e. the quantum limit SC correlation function), and setting the tuning parameters to infinity results in a classical limit description of the system similar to LSC-IVR. MQC-IVR thus provides a uniform framework for SC simulations with mode-specific quantization and no uncontrolled approximations to the forces between quantum and classical subsystems, as are made in standard multi-physics approaches.Caro and Salcedo (1999)

A challenge with MQC-IVR is, however, determining the optimal set of tuning parameters to minimize both computational cost and loss of accuracy. One approach to this challenge, and the subject of this study, is to analytically evaluate the general MQC-IVR limit where the tuning parameters associated with the quantum subsystem go to zero, and the tuning parameters associated with the classical subsystem go to infinity. The result is an analytical mixed quantum-classical (AMQC-IVR) expression for the SC-IVR time correlation function that offers reduced computational effort, and circumvents the need to find optimal values of the tuning parameters used in MQC-IVR.

In this study we use three multidimensional model systems to demonstrate that AMQC-IVR accurately describes quantum dynamical features of systems over a wide range of coupling strengths between quantum and classical subsystems. We also introduce a separable prefactor (SP) approximation to the AMQC-IVR prefactor that further reduces computational cost. We show that the SP approximation is increasingly accurate as the coupling between quantum and classical subsystems decreases, and increasingly efficient when the classical subsystem is larger than the quantum subsystem. Finally, we show that AMQC-IVR has the potential to be systematically improved, and to be amenable to a variety of existing representations and approximations in the SC literature.Gelavert et al. (2000); Di Liberto and Ceotto (2016)

This paper is organized as follows. In Sec. II we briefly review MQC-IVR theory and provide an overview of the derivation of AMQC-IVR. In Sec. III we describe the model systems and in Sec. IV we provide simulation details. Results are discussed in Sec. V and conclusions are drawn in Sec. VI.

II Theory

II.1 MQC-IVR

Throughout this manuscript we use atomic units and take $\hbar=1$. The general MQC-IVR correlation functionChurch, Antipov, and Ananth (2017); Church et al. (2018) is given by

[TABLE]

where $N$ is the dimensionality of the entire system, and we use the double-forward formulation. Church, Antipov, and Ananth (2017) The $2N$-dimensional phase space vectors of the unprimed and primed trajectories at time $t$ are defined as

[TABLE]

respectively. We also take $S_{t}\left(\mathbf{z}\right)$ to be the classical action of a trajectory originating at point $\mathbf{z}$, and $\bm{\Delta}_{0}=\mathbf{z}_{0}^{\prime}-\mathbf{z}_{0}$ is the phase space displacement between pairs of forward trajectories at $t=0$. The position-space wavefunction of the coherent state $\ket{\mathbf{z}_{t}}$ at time $t$ is given by

[TABLE]

where $\bm{\gamma}_{t}$ is a diagonal $N\times N$ matrix that determines the width of the wavepacket. We also represent the coherent state matrix element of a given quantum mechanical operator $\hat{\Omega}$ as

[TABLE]

A detailed form of the MQC-IVR prefactor $D_{t}\left(\mathbf{z}_{0},\mathbf{z}_{0}^{\prime};\mathbf{c}\right)$ is provided in Appendix A. The diagonal $2N\times 2N$ matrix of tuning parameters $\mathbf{c}$ is given by

[TABLE]

and $\mathbb{O}$ is the null matrix. The elements of the diagonal $N\times N$ matrices $\mathbf{c}_{\mathbf{p}}$ and $\mathbf{c}_{\mathbf{q}}$ determine the extent of separation (in momentum and position space, respectively) between trajectory pairs at time $t=0$ and, therefore, determine the extent of phase cancellation in the integrand of Eq. 1. In the limit that all the elements of $\mathbf{c}$ approach zero, the effect of MFF is removed and Eq. 1 reduces to DHK-IVR,

[TABLE]

which we refer to as the quantum limit SC-IVR correlation function. In Eq. 7, $C_{t}\left(\mathbf{z}_{0}\right)$ is the Herman-Kluk prefactor for a trajectory beginning at point $\mathbf{z}_{0}$,

[TABLE]

and elements of the monodromy matrix $\mathbf{M}$ are defined by $\mathbf{M}_{\alpha\beta}=\frac{\partial\bm{\alpha}_{t}}{\partial\bm{\beta}_{0}}$ with $\left(\bm{\alpha},\bm{\beta}\right)\in\left(\mathbf{p},\mathbf{q}\right)$. In the limit that all the elements of $\mathbf{c}$ approach infinity, the unprimed and primed trajectories are constrained to be identical, resulting in complete phase cancelation in the MQC-IVR integrand. In this limit, MQC-IVR is identical to a classical limit SC-IVR correlation function similar to LSC-IVR,

[TABLE]

which we refer to as Husimi-IVR. It is clear that the values of $\mathbf{c}$ in MQC-IVR control the phase contributions to the integrand from each dof and, therefore, the level of SC theory used to describe each dof.

II.2 AMQC-IVR

Here we outline the AMQC-IVR derivation and provide details in Appendix B. First, as a matter of bookkeeping, we consider a general system with $F$ quantum dofs and $N-F$ classical dofs. We also order the elements of the position and momentum vectors ($\mathbf{q}_{t}$, $\mathbf{q}^{\prime}_{t}$, $\mathbf{p}_{t}$, and $\mathbf{p}^{\prime}_{t}$) such that the $F$ elements of the quantum subsystem are listed before the $N-F$ elements of the classical subsystem.

The $\mathbf{c}$-dependence of the MQC-IVR integrand in Eq. 1 can be written as

[TABLE]

We now evaluate Eq. 10 in the limit that the elements of $\mathbf{c}$ associated with the quantum subsystem go to zero, and the limit that the elements of $\mathbf{c}$ associated with the classical subsystem go to infinity. In order to distinguish the phase space variables of the two subsystems, we introduce the following $2F$-dimensional vectors to represent the initial conditions of the quantum subsystem,

[TABLE]

and the following $2(N-F)$-dimensional vectors to represent the initial conditions of the classical subsystem,

[TABLE]

Using the following $\delta$-function identity,

[TABLE]

we obtain

[TABLE]

Note that the product on the right-hand side of Eq. 16 is over all $2\times(N-F)$ elements of the classical subsystem in $\mathbf{z}_{C}$ and $\mathbf{z}_{C}^{\prime}$. In addition, we note that the prefactor can be written as a product of a pair of Herman-Kluk prefactors for the primed and unprimed trajectories and $\Lambda_{t}\left(\mathbf{z}_{0},\mathbf{z}_{0}^{\prime}\right)$, an additional term that accounts for coupling between the quantum and classical subsystems defined in Appendix B. Substituting the limit of Eq. 16 into Eq. 1, and then evaluating the integrals over the primed initial conditions of the classical subsystem, i.e. over $d\mathbf{z}_{C}^{\prime}$, we obtain the AMQC-IVR correlation function,

[TABLE]

The $2N$-dimensional phase space vector $\bar{\mathbf{z}}_{0}$ replaces the classical components of $\mathbf{z}_{0}^{\prime}$ with the classical components of $\mathbf{z}_{0}$.

In obtaining Eq. 17 from Eq. 1 we have constrained the initial conditions of the classical subsystem in the unprimed and primed trajectories to be identical, and, consequently, reduced the dimensionality of the phase space integral by $2\times(N-F)$, twice the dimensionality of the classical subsystem. The additional prefactor $\Lambda_{t}$ appears as a result of this constraint. The reduced dimensionality of the integral and the similarity of initial conditions in each trajectory pair also results in significant phase cancellation and an acceleration in convergence, while the initial non-zero displacements in the quantum subsystem contribute phase information essential to describing quantum coherence effects accurately.

Like MQC-IVR, AMQC-IVR is an approximate representation of the quantum mechanical time correlation function that offers mode-specific quantization in a dynamically uniform framework. While both methods offer a significant improvement over DHK-IVR’s sign problem, AMQC-IVR does not require any tuning parameters, making it a more efficient and theoretically satisfying implementation of the mixed SC limit. Furthermore, as discussed in the following section, the AMQC-IVR prefactor is a good starting point for a variety of additional approximations that can improve computational efficiency even further.

II.3 Separable Prefactor Approximation

In the simplest case where the quantum and classical subsystems are not coupled, AMQC-IVR can be simplified significantly. First, note that the phase terms and prefactors in the MQC-IVR and DHK-IVR integrands can be separated into the product of a quantum component and a classical component,

[TABLE]

Matrices $\mathbf{c}_{Q}$ and $\mathbf{c}_{C}$ represent the quantum and classical blocks of matrix $\mathbf{c}$. Given the limiting behavior of MQC-IVR as described in Sec. II, and given the separability of Eq. 18, it is straightforward to show that the AMQC-IVR expression for non-interacting quantum-classical subsystems is given by

[TABLE]

Note that the phase terms and prefactors in Eq. 20 only depend on the quantum subsystem. A comparison of Eq. 20 and Eq. 17 with the use of Eq. 19 shows that, when quantum and classical subsystems do not interact, the AMQC-IVR prefactor reduces to unity,333Note that, since the unprimed and primed trajectories of the classical subsystem are identical (when uncoupled from the quantum subsystem), the product of phase terms and prefactors associated with the classical subsystem is also equal to unity: $C_{t}\left(\mathbf{z}_{C}\right)C_{t}^{*}\left(\mathbf{z}_{C}\right)e^{i\left[S_{t}\left(\mathbf{z}_{C}\right)-S_{t}\left(\mathbf{z}_{C}\right)\right]}=1$.

[TABLE]

We now use Eq. 20 and Eq. 21 to motivate an efficient approximation to AMQC-IVR when quantum and classical subsystems are weakly coupled,

[TABLE]

We refer to Eq. 22 as the SP approximation, and it is increasingly valid as the coupling between quantum and classical subsystems approaches zero. The advantage of the SP approximation is that it significantly reduces computational expense, particularly when the classical subsystem is larger than the quantum subsystem, i.e. when $N>>F$. This is because the SP approximation contains only two $F\times F$ Herman-Kluk prefactors, rather than two $N\times N$ Herman-Kluk prefactors and the additional $4N\times 4N$ prefactor $\Lambda_{t}\left(\mathbf{z}_{0},\mathbf{z}_{Q}^{\prime}\right)$ in Eq. 17.

There are several other approximations that could be made in addition to the SP approximation. For example, one could choose to evolve only a subset of the monodromy matrix elements, those associated with the quantum subsystem, albeit approximately, and reduce computation time even further. Furthermore, since the SP approximation contains only Herman-Kluk prefactors, which have been extensively studied, it is amenable to a variety of other representations and approximationsGelavert et al. (2000); Di Liberto and Ceotto (2016) that may reduce computational expense even further. In this study, we numerically explore the applicability of the SP approximation, reserving other possibilities for future work.

III Model Systems

Model 1 is a 1D anharmonic oscillator coupled to a heavy harmonic ‘bath’ mode. The Hamiltonian is given by

[TABLE]

where $m_{i}$ and $\omega_{i}$ are the mass and frequency of the $i^{th}$ mode and $k$ is the bilinear coupling parameter. The initial state of the system is a product of coherent states, the position-space wavefunction of which is

[TABLE]

with $\mathcal{N}$ for normalization.

Model 2 contains the anharmonic mode of model 1 but now coupled to a bath of $N-1$ harmonic oscillators. The Hamiltonian is given by

[TABLE]

The initial state of the full system is a product of $N$ coherent states, as in Eq. 25.

Model 3 is the widely used 1D symmetric double-well potential coupled to a thermal bath of $N-1$ harmonic oscillators. The Hamiltonian is

[TABLE]

with

[TABLE]

In models 2 and 3, an Ohmic spectral density with an exponential cutoff is used for the bath,

[TABLE]

In model 2 we take $\omega_{c}=\omega_{1}$ and in model 3 we take $\omega_{c}=\omega_{b}$.

IV Simulation Details

With model 1 we take $\hat{B}=\hat{x}_{1}$ to compute the average position of the anharmonic mode $\langle\hat{x}_{1}\rangle_{t}$ as a function of time. We take $m_{1}=1.0$, $m_{2}=25.0$, $\omega_{1}=\sqrt{2}$, $\omega_{2}=1/3$, and we vary the bilinear coupling $k$ between $0.5$ and $2.0$, all in atomic units. Operator $\hat{A}$ is the projection operator corresponding to the initial state,

[TABLE]

The initial coherent states are centered at $q_{i_{1}}=q_{i_{2}}=1.0$ and $p_{i_{1}}=p_{i_{2}}=0.0$ with width parameters $\gamma_{j}=m_{j}\omega_{j}$ $\left(j=1,2\right)$. The initial conditions of the $j^{th}$ quantum dof are sampled from

[TABLE]

and the initial conditions of the $j^{th}$ classical dof are sampled from

[TABLE]

with $\mathcal{N}$ for normalization. The coherent state matrix element of $\hat{x}_{1}$ is

[TABLE]

We found that approximately $10^{5}$ trajectory pairs are needed to converge AMQC-IVR calculations with model 1, whereas DHK-IVR required upwards of $10^{8}$. We also use a time step of $0.05$.

With model 2 we compute $\langle\hat{x}_{1}\rangle_{t}$ in the presence of $12$ bath modes as a function of time, and with different values of the coupling $\eta$ between the system and the bath. The initial coherent states are centered at $q_{i_{1}}=1.0$, $p_{i_{1}}=0.0$, and $q_{i_{j}}=p_{i_{j}}=0.0$ $\forall$ $j\in\left[2,N=13\right]$. And the coherent state width parameter of the $j^{th}$ dof is given by $\gamma_{j}=m_{j}\omega_{j}$. The initial conditions of the quantum dof are sampled from Eq. 31, and the initial conditions of the classical dofs are sampled from Eq. 32. We found that AMQC-IVR requires approximately $10^{6}$ trajectory pairs for convergence. We also use a time step of $0.025$.

With model 3 we compute the thermal transmission coefficient $\kappa\left(T\right)$ of the symmetric double-well by means of a flux-side correlation function,Miller (1993)

[TABLE]

where $k_{\text{cl}}^{\text{TST}}(T)$ is the classical transition state theory result, $Q_{r}(T)$ is the partition function in the reactant well, and $C_{fs}(t)$ is the flux-side correlation function characterized by the following operators,

[TABLE]

Eq. 36 contains the flux operator $\hat{F}=i\left[\hat{H},\hat{h}\right]$ and constant $\beta=1/kT$ with temperature $T=300K$ and Boltzmann constant $k$. In Eq. 37 $\hat{h}$ is the unit step function specifying the dividing surface. Numerical parameters for the Hamiltonian of model 3 can be found in Ref. [22]. In order to evaluate the coherent state matrix element of $\hat{A}$ in Eq. 36 we make a normal-mode approximation at the transition stateWang, Thoss, and Miller (2000) so that the Hamiltonian is approximately separable:

[TABLE]

Frequencies $\lambda^{\ddagger}$ and $\lambda_{j}$ are the imaginary and real normal-mode frequencies at the transition state, respectively. Under this approximation, the coherent state matrix element of $\hat{A}$ in Eq. 36 is given by

[TABLE]

with $u_{1}=\beta\left|\lambda^{\ddagger}\right|/2$, $u_{j}=\beta\lambda_{j}/2$ $\forall$ $j$ $=$ $2\dots N$, and we choose $\gamma_{1}=m_{1}|\lambda^{\ddagger}|\cot u^{\ddagger}$ as well as $\gamma_{j}=m_{j}\lambda_{j}$ $\forall$ $j$ $=$ $2\dots N$. With AMQC-IVR, however, only the diagonal elements of Eq. 43 are needed, i.e.

[TABLE]

Note that the products on the right-hand sides of Eq. 41 and Eq. 44 are over all components of the classical subsystem. Initial conditions for the quantum and classical dofs are sampled from the following distributions,

[TABLE]

respectively. Constant $\mathcal{N}$ is for normalization. The coherent state matrix element of $\hat{h}$ at time $t$ is given by

[TABLE]

$\text{Erf}\left(x\right)$ is the error function of $x$. We found that approximately $10^{6}$ trajectory pairs were required for convergence. We also use time step of $10.0$. The SC Corr-Code Package, an open-source program developed in-house, was used to run all the simulations in this study.cor

V Results and Discussion

The position expectation value of the anharmonic mode in model 1 is plotted as a function of time in Fig. 1 with three different coupling strengths. The rapid oscillations in the exact quantum results (dashed) at long times are a result of nuclear coherence, a feature that is clearly absent in the classical limit Husimi-IVR results (black). The AMQC-IVR results (blue), in which the anharmonic mode is treated in the quantum limit and the harmonic mode is treated in the classical limit, is consistently accurate (see Table 1) with only slight damping in the coherences at long times when the quantum and classical modes are strongly coupled.

Also plotted in Fig. 1 are the AMQC-IVR results obtained with the SP approximation (red) of Eq. 22. In Fig. 1(a) and Fig. 1(b) these results are nearly identical to the AMQC-IVR results without approximations to the prefactor. There is some additional loss of long-time coherences in the SP approximation in Fig. 1(c), and the reported error increases (see Table 2). We note that the SP results with this low-dimensional model serve only as a proof of principle, offering limited reduction in computational cost.

Fig. 1(c) also plots $\langle\hat{x}_{1}\rangle_{t}$ as computed with AMQC-IVR when both modes are treated in the quantum limit. Since all dofs are quantized here, this is equivalent to using DHK-IVR of Eq. 7, but the AMQC-IVR prefactors were used during computation. The result is nearly identical to the exact quantum result at all times, showing that, by increasing the size of the quantum subsystem, AMQC-IVR results can be systematically improved toward the DHK-IVR limit.

The average position of the anharmonic mode of model 2 is plotted in Fig. 2 as a function of time. Each panel uses a different reduced coupling strength $\eta/m\omega_{1}$ between the anharmonic dof and the harmonic bath. And each AMQC-IVR result was obtained by quantizing the anharmonic mode and treating the harmonic bath in the classical limit.

In Fig. 2(a), where the coupling between the system and bath is very weak, AMQC-IVR (blue) is nearly identical to the exact result of the uncoupled anharmonic oscillator (dashed), but with a slight overestimation of the amplitudes at later times. The SP approximation (red) in Fig. 2(a) is nearly identical the exact uncoupled result as well, and computational time was reduced by a factor of four. The classical limit result obtained with Husimi-IVR (black) in Fig. 2(a) appears to fail at all but very short times. A comparison of these results suggest that AMQC-IVR and the SP approximation are accurately capturing nuclear coherence effects in the position correlation function. This is encouraging since a full DHK-IVR treatment of such a highly multidimensional system is not possible without a remedy to the sign problem.

In Fig. 2(b), where the coupling between the system and bath is $100$ times stronger than in Fig. 2(a), the AMQC-IVR and SP results again resemble the oscillatory structure of the uncoupled result, albeit with damped amplitudes around $t=65$ due to the stronger influence from bath dofs. Once again, the Husimi-IVR result in Fig. 2(b) does not show the oscillatory behavior of the AMQC-IVR with and without the SP approximation.

Fig. 2(c), where the coupling is $10^{4}$ times stronger than in Fig. 2(a), the AMQC-IVR and Husimi-IVR results are very similar, though the amplitudes in the AMQC-IVR result are slightly larger from about $t=20$ and beyond. Since the coupling between the system and bath is strong here, one would expect that the coherence effects associated with the motion of the system would be significantly damped due the strong influence from the bath. It is therefore reasonable to see that AMQC-IVR and Husimi-IVR give similar results in the strong coupling limit. The oscillatory structure of the SP approximation in Fig. 2(c) resembles that of the AMQC-IVR and Husimi-IVR results, but with some moderate damping of the amplitudes, particularly at longer times. However, since the SP approximation is most valid when the coupling between the system and bath is weak, it is intuitive that the SP approximation is less reliable than AMQC-IVR in Fig. 2(c).

The thermal transmission coefficient of the symmetric double well of model 3 is plotted in Fig. 3 as a function of the reduced coupling strength $\eta/m\omega_{b}$ (see Appendix C for the tabulated data). AMQC-IVR results (red) were obtained by quantizing the double-well and treating each bath dof in the classical limit.

AMQC-IVR clearly agrees well with the exact quantum results (black) in the weak coupling limit, and captures the turnover region around $\eta/m\omega_{b}=0.5$. In the stronger coupling limit, however, AMQC-IVR begins to increasingly underestimate the transmission coefficient, though the trend is qualitatively correct. It is reasonable to expect AMQC-IVR to fail in the strong coupling limit, since much of the phase information from bath dofs is removed.

Fig. 4 plots the thermal transmission coefficient of model 3 in the weak coupling regime, $\eta/m\omega_{b}=0.05$, as a function of time with AMQC-IVR (black) and the SP approximation (red).

Both results show an oscillatory structure in the transmission coefficient, corresponding to the transfer of population between the two wells. The amplitude of the second peak ($t=6000$) in the SP result is slightly damped relative to the AMQC-IVR result, but, appealingly, both results clearly converge to the same long-time limit, where the rate is determined. Furthermore, computation with the SP approximation was $4$ times faster than without, and when the number of bath modes was doubled to $24$ (which does not affect the rate), computation was $8$ times faster than without. Note, however, that while the $\%$ error between AMQC-IVR and the SP approximation is small when $\eta/m\omega_{b}=0.05$, the SP approximation increasingly underestimates the rate for larger coupling strengths (see Table 3).

Overall, our analysis of model 3 shows that AMQC-IVR and the SP approximation can be very reliable in computing chemical reaction rates in the condensed phase, particularly when the coupling between quantum and classical subsystems is weak.

VI Conclusions

An analytical parameter-free mixed quantum-classical limit of DHK-IVR has been derived such that some modes of the system can be treated in the quantum and others in the classical limit of SC-IVR theory. AMQC-IVR shares the advantages of MQC-IVR as a methodology that both mitigates the SC-IVR sign problem and offers mode-specific quantization in a dynamically uniform framework. But AMQC-IVR offers a reduction in dimensionality of the phase space integral and it removes the task of having to choose an optimal set of tuning parameters. We have tested AMQC-IVR on three system-bath models and showcased its numerical accuracy across coupling regimes. We have also shown that the SP approximation to AMQC-IVR can be very accurate when the coupling between quantum and classical dofs is weak, and very efficient when the quantum subsystem is smaller than the classical subsystem.

Acknowledgements

The authors thank Prof. Makri for providing the exact quantum results for the transmission coefficient used in Fig. 3. This work was funded by an NSF CAREER Grant No. CHE 1555205 and the Research Corporation for Science Advancement through a Cottrell Scholar Award.

Appendix A MQC-IVR Prefactor

The MQC-IVRChurch, Antipov, and Ananth (2017) prefactor is given by

[TABLE]

with

[TABLE]

Appendix B Derivation of AMQC-IVR

Here we derive AMQC-IVR in detail. According to a previous studyChurch, Antipov, and Ananth (2017) one can rewrite Eq. 16 as,

[TABLE]

$\mathbf{K}^{T}$ is a complex $4N\times 4N$ matrix,

[TABLE]

with constant $2N\times 2N$ matrices

[TABLE]

$\mathbf{M}$ and $\mathbf{M}^{\prime}$ are the full $2N\times 2N$ monodromy matrices for trajectories beginning at $\mathbf{z}_{0}$ and $\mathbf{z}_{0}^{\prime}$, respectively. We also have the $4N\times 4N$ diagonal matrix of tuning parameters,

[TABLE]

with $\mathbf{c}$ defined in Eq. 6. We also have the $4N\times 4N$ matrix $\mathbf{J}$,

[TABLE]

The procedure from here is to expand the determinant on the right-hand side of Eq. 57 and evaluate the limit term-by-term. Using the definition of the determinant of a general $4N\times 4N$ matrix $\mathbf{\Xi}$,

[TABLE]

where $\epsilon$ is the Levi-Civita symbol, the expansion of the determinant in Eq. 57 gives

[TABLE]

Given that $\tilde{\mathbf{c}}$ is diagonal, and given that the lower blocks of $\tilde{\mathbf{c}}$ contain only zeros, we have

[TABLE]

When the leading product on the right-hand side of Eq. 65 is carried through the summation, the only terms that will survive the limit are those that are independent of the elements of $\mathbf{c}$. We then have,

[TABLE]

Note that the right-most product on the right-hand side of Eq. 66 is over all $m$ $\in\left[1,F\right]$, all $m$ $\in\left[N+1,N+F\right]$, and all $m$ $\in\left[2N+1,4N\right]$. The right-hand side Eq. 66 can be viewed as the determinant of matrix $\mathbf{K}^{T}$ plus an additional matrix $\mathbf{\Sigma}$,

[TABLE]

The addition of matrix $\mathbf{\Sigma}$ to matrix $\mathbf{K}^{T}$ effectively replaces certain rows of $\mathbf{K}^{T}$ (specifically, rows associated with the classical subsystem) with the corresponding rows of matrix $i\mathbf{J}$. Therefore,

[TABLE]

Since matrix $\mathbf{K}^{T}$ is invertible, we can equivalently write Eq. 67 as

[TABLE]

And since the leading term on the right-hand side of Eq. 70 is proportional to the product of HK prefactors,

[TABLE]

the AMQC-IVR prefactor is given by,

[TABLE]

Note that, for the model systems studied here, however, Eq. 67 was used for all computations. The results of Eq. 70, Eq. 71, and Eq. 72 are used here to motivate the SP approximation.

Appendix C Tabulated Thermal Transmission Coefficients