Thermal quantum time-correlation functions from classical-like dynamics

TL;DR

This paper reviews recent advances in classical-like dynamics methods such as CMD, RPMD, and TRPMD for calculating thermal quantum time-correlation functions, emphasizing their derivation from Matsubara dynamics and their relation to quantum transition-state theory.

Contribution

It provides a unified derivation of these methods from Matsubara dynamics and clarifies their connection to quantum transition-state theory.

Findings

Methods like CMD, RPMD, and TRPMD are derived from Matsubara dynamics.

Matsubara-TST is shown to be equivalent to RPMD-TST.

The paper identifies future directions for improving quantum dynamics simulations.

Abstract

Thermal quantum time-correlation functions are of fundamental importance in quantum dynamics, allowing experimentally-measurable properties such as reaction rates, diffusion constants and vibrational spectra to be computed from first principles. Since the exact quantum solution scales exponentially with system size, there has been considerable effort in formulating reliable linear-scaling methods involving exact quantum statistics and approximate quantum dynamics modelled with classical-like trajectories. Here we review recent progress in the field with the development of methods including Centroid Molecular Dynamics (CMD), Ring Polymer Molecular Dynamics (RPMD) and Thermostatted RPMD (TRPMD). We show how these methods have recently been obtained from `Matsubara dynamics', a form of semiclassical dynamics which conserves the quantum Boltzmann distribution. We also rederive t->0+ quantum transition-state theory (QTST) in the Matsubara dynamics formalism showing that Matsubara-TST, like RPMD-TST, is equivalent to QTST. We end by surveying areas for future progress.

Tables

Table 1. Table 1: Summary of the properties of LSC-IVR, CMD, RPMD and TRPMD. T c subscript 𝑇 𝑐 T_{c} is the crossover temperature discussed in appendix  D .

	LSC-IVR	CMD	RPMD	TRPMD
Approximation	Discard $𝒪\left({\hslash }^2\right)$	Mean field	Discard $i{ℒ}_{\Im }^{\left[M\right]}$	Replace $i{ℒ}_{\Im }^{\left[M\right]}$ with ${𝒜}_{\mathrm{wn}}^{\left[M\right]†}$
Conserves distribution and detailed balance?	No	Yes	Yes	Yes
Centroid force	N/A	Mean field	Matsubara force	Matsubara force
Reaction rates	Problems beneath $T_c$[50]	Inaccurate beneath $T_c$ [77, 56]	Good [4, 5]	Friction slows rates [23]
Spectra	Good[49]	Curvature problem [34, 72]	Spurious resonances [34, 72]	Good [70]
Diffusion	ZPE leakage [49]	Good [70, 78]	Good [49]	Good [79]
Nonlinear operators	Good if ZPE not problematic [6]	Fails even at $t=0$ [60]	Breakdown from incorrect frequencies [60]	Breakdown from damping [59]
Advised usage	Nonlinear operators	Rates above $T_c$, diffusion	Rates, diffusion	Spectra, diffusion

Equations

$$H(p,q)=\frac{p^2}{2m}+V(q).$$

$$G_{AB}(t)=\frac{1}{2\pi \hslash }\int dp\int dq{e}^{-\beta H(p,q)}A(p,q)B(p_t,q_t)$$

$$G_{AB}(t)=\frac{1}{2\pi \hslash }\int dp\int dq{e}^{-\beta H(p,q)}A(p,q)B(p,q,t)$$

$$\frac{d}{dt}B(p,q,t)=$$

$$=$$

$${\left(\frac{\partial B(p,q,t)}{\partial t}\right)}_{p,q}=0,$$

$$\mathcal{L}=\frac{p}{m}\frac{\partial }{\partial q}-\frac{\partial V(q)}{\partial q}\frac{\partial }{\partial p}$$

$$\frac{d}{dt}B(p,q,t)=\mathcal{L}B(p,q,t)$$

$$G_{AB}(t)=\frac{1}{2\pi \hslash }\int dp\int dq{e}^{-\beta H(p,q)}A(p,q)e^{\mathcal{L}t}B(p,q,0).$$

$$\frac{d}{dt}B(p_t,q_t)=\frac{\partial B(p_t,q_t)}{\partial q_t}\frac{d{q}_t}{dt}+\frac{\partial B(p_t,q_t)}{\partial p_t}\frac{d{p}_t}{dt}$$

$$\mathcal{L}B(p_t,q_t)=\frac{\partial B(p_t,q_t)}{\partial q_t}\mathcal{L}{q}_t+\frac{\partial B(p_t,q_t)}{\partial p_t}\mathcal{L}{p}_t.$$

$$\frac{d{q}_t}{dt}=\mathcal{L}{q}_t,\frac{d{p}_t}{dt}=\mathcal{L}{p}_t$$

$$\frac{d}{dt}{G}_{AB}(t)=-\frac{1}{2\pi \hslash }\int dp\int dq{e}^{-\beta H(p,q)}A(p,q)\overleftarrow{\mathcal{L}}B(p,q,t)$$

$$\frac{d}{dt}{G}_{AB}(t)=-\frac{1}{2\pi \hslash }\int dp\int dq{e}^{-\beta H(p,q)}B(p,q,t)\mathcal{L}A(p,q).$$

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

$$=$$

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

$$\hat{H}=\frac{{\hat{p}}^2}{2m}+V(\hat{q}).$$

$$c_{AB}(t)=\mathrm{Tr}\left[e^{-\beta \hat{H}}\hat{A}{e}^{i\hat{H}t/\hslash }\hat{B}{e}^{-i\hat{H}t/\hslash }\right]$$

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

$$\times ⟨y|e^{i\hat{H}t/\hslash }|z⟩B(z)⟨z|e^{-i\hat{H}t/\hslash }|x⟩$$

$$c_{AB}(-t{)}^{\ast }=c_{BA}(t).$$

$${\bar{c}}_{AB}(t)=:$$

$$=$$

$${\bar{c}}_{AB}(t)=$$

$$\times ⟨q+\Delta /2|e^{i\hat{H}t/\hslash }|z⟩B(z)⟨z|e^{-i\hat{H}t/\hslash }|q-\Delta /2⟩.$$

$${\tilde{c}}_{AB}(t)=\frac{1}{\beta }{\int }_0^{\beta }d\lambda \mathrm{Tr}[e^{-(\beta -\lambda )\hat{H}}\hat{A}{e}^{-\lambda \hat{H}}{e}^{i\hat{H}t/\hslash }\hat{B}{e}^{-i\hat{H}t/\hslash }]$$

$${\tilde{c}}_{AB}(t)={\tilde{c}}_{AB}(t{)}^{\ast }$$

$${\tilde{c}}_{AB}(-t)={\tilde{c}}_{BA}(t)$$

$$C_{AB}^{[N]}(t)=\int d\mathbf{q}\int d\mathbf{\Delta }$$

$$\times \prod _{i=0}^{N-1}⟨q_{i-1}-{\Delta }_{i-1}/2|\frac{1}{2}(\hat{A}{e}^{-{\beta }_N\hat{H}}+e^{-{\beta }_N\hat{H}}\hat{A})|q_i+{\Delta }_i/2⟩$$

$$\times ⟨q_i+{\Delta }_i/2|e^{i\hat{H}t/\hslash }\hat{B}{e}^{-i\hat{H}t/\hslash }|q_i-{\Delta }_i/2⟩$$

$$\hat{A}=\frac{1}{N}\sum _{k=0}^{N-1}{\hat{A}}_k$$

$$\underset{N\to \infty }{\lim }{C}_{AB}^{[N]}(t)={\tilde{c}}_{AB}(t).$$

Full text

1 Introduction

Quantum thermal time-correlation functions [1, 2] are routinely used to calculate reaction rates, spectra and diffusion constants amongst many other physically observable quantities, and provide a useful bridge between the algebra of quantum mechanics and experimental measurement. In general they can only be computed exactly for very small or model systems, and there is consequently a need for reliable approximate computation with classical-like scaling (i.e. linear scaling w.r.t. the number of dimensions of the system). The purpose of this New View article is to review the origins of a number of these methods; namely the approximations they make to the exact quantum evolution and the conditions under which they are likely to be valid. This should allow a theoretician to discern for themselves the optimal method for a given problem.

This article is designed to provide an overview of the field with references for further reading and is not intended to be exhaustive. Applications of many of the methods discussed here have already been extensively reviewed, including centroid molecular dynamics (CMD) [3], ring polymer molecular dynamics (RPMD) [4], RPMD rate theory [5] and the linearized semiclassical initial-value representation (LSC-IVR) [6]. Consequently, applications of these methods are only mentioned when pertinent.

We also rederive quantum transition-state theory (QTST) in the Matsubara formalism, showing that Matsubara TST is identical to QTST provided that the dividing surface is only a function of the Matsubara modes, and which in turn is identical to RPMD-TST when the dividing surface is invariant to cyclic permutation in imaginary time. For reviews on rate theory more generally, see Refs. [7, 8, 9, 10].

There exist many other methods to simulate quantum dynamics which are not covered here, including exact quantum methods such as multi-configuration time-dependent Hartree (MCTDH) [11], matrix-based methods [12], and path-integrals [13]. Other approaches include gaussian wavepacket propagation [14], semiclassical dynamics [15, 16] and mixed quantum-classical dynamics [17, 18, 19].

For most of the article we assume that dynamics is on a single Born-Oppenheimer potential energy surface that is known and differentiable (either of a model form, fitted to some set of parameters, or from ab initio electronic structure theory); the computation of accurate potential energy surfaces is a discipline in itself. We touch upon extensions to non-adiabatic dynamics towards the end. We generally assume that the systems being described are in thermal equilibrium; application to non-equilibrium systems is an interesting area of present research [20].

The article is structured as follows. In section 2 we review classical and quantum thermal time-correlation functions, the Wigner transform and the Moyal series. Section 3 touches upon LSC-IVR, and section 4 provides the derivation of Matsubara dynamics. Section 5 covers approximations to Matsubara dynamics such as CMD, RPMD and TRPMD, and section 6 gives an alternative derivation of QTST in the Moyal/Matsubara formalism. Section 7 presents directions for future research and section 8 concludes.

2 Thermal time-correlation functions

Here we briefly present background theory sufficient to follow the remainder of the article; further detail is available in standard texts[21, 2, 1].

2.1 Classical

For simplicity we consider a one-dimensional system, extension to further dimensions being straightforward[2], with position $q$ and momentum $p$ and a classical Hamiltonian

[TABLE]

The thermal correlation function between observables $A$ and $B$ at inverse temperature $\beta\equiv 1/k_{\rm B}T$ (where $k_{\rm B}$ is the Boltzmann constant) is generally written as

[TABLE]

where $p$ and $q$ are sampled at zero time and $(p_{t},q_{t})\equiv(p_{t}(p,q,t),q_{t}(p,q,t))$ are the solutions to a classical trajectory for length $t$ starting at $(p,q)$ at time $t=0$. The correlation function can equivalently be given as

[TABLE]

where $B(p,q,t)$ corresponds to an initial phase-space distribution $(p,q)$ propagated for time $t$. Formally, one can obtain the dynamical equations of motion by differentiating Eq. (3) w.r.t. time to obtain

[TABLE]

where we have applied Newton’s first and second law to obtain Eq. (5). Strictly speaking, we are also assuming that the observables themselves are not explicit functions of time, i.e.

[TABLE]

and likewise for $A$, which is the case for all correlation functions considered in this article. Equation (5) allows us to define a classical Liouvillian111Following the convention of Zwanzig [1] we define the Liouvillian without a prefactor of $i$.

[TABLE]

such that

[TABLE]

which has a formal solution $B(p,q,t)=e^{\mathcal{L}t}B(p,q,0)$ and therefore

[TABLE]

To see how Eq. (2) is equivalent to Eq. (3) we differentiate Eq. (2) w.r.t. $t$, obtaining

[TABLE]

but if Eq. (2) is a solution to Eq. (3) then by Eq. (8), the LHS of Eq. (10) must be equal to the action of the Liouvillian on $B(p_{t},q_{t})$, which is

[TABLE]

Comparing Eq. (10) and Eq. (11) gives

[TABLE]

which have formal solutions $q_{t}=e^{\mathcal{L}t}q,\ p_{t}=e^{\mathcal{L}t}p$. This means that instead of propagating a phase space density in $B(p,q,t)$, one can simply propagate individual positions and momenta to find $(p_{t},q_{t})$ and insert into the function $B(p_{t},q_{t})$, which is computationally easier. However, if $\mathcal{L}$ contains higher derivatives in $p$ and/or $q$ (as is the case in exact quantum evolution and stochastic dynamics) then this convenient property no longer holds.

If $B=H$, then from Eq. (7) $\mathcal{L}H=0$, meaning that classical dynamics conserves the classical Hamiltonian, as to be expected. It follows that $\mathcal{L}e^{-\beta H(p,q)}=0$ and the classical dynamics conserves the classical Boltzmann distribution.

If we differentiate Eq. (3) w.r.t. $t$, apply Eq. (7) use integration by parts on the derivatives in $p$ and $q$ we obtain

[TABLE]

where $\overleftarrow{\mathcal{L}}$ is ‘acting backwards’ onto $e^{-\beta H(p,q)}A(p,q)$, but using the product rule and that $\mathcal{L}e^{-\beta H(p,q)}=0$, this gives

[TABLE]

Integration of this, noting that $B(p,q,0)=B(p,q)$ gives

[TABLE]

which is detailed balance. Note that this is a stronger condition than time reversal symmetry, which only implies [from Eq. (13)] that

[TABLE]

where the distribution has to be propagated too. In general, if the dynamics conserves the distribution then the correlation function will observe detailed balance (strictly speaking, for stochastic systems this is a necessary but not sufficient requirement [22, 23]).

2.2 Quantum

Similar to the classical case, we consider a one-dimensional system with mass $m$, co-ordinate $q$ with conjugate momentum $p$ and quantum Hamiltonian

[TABLE]

In this section we introduce a variety of quantum time-correlation functions and briefly discuss their properties, particularly concerning the ease with which they may be approximated by classical-like dynamics.

2.2.1 Conventional time-correlation function

The conventional quantum time-correlation function is given by[2, 24]

[TABLE]

such that $c_{AB}(0)={\rm Tr}[e^{-\beta\hat{H}}\hat{A}\hat{B}]$, giving the thermal average of $\hat{A}$ and $\hat{B}$. Since $[e^{-i\hat{H}t/\hbar},\hat{H}]=0$, $c_{AH}(t)=c_{AH}(0)$ and the quantum dynamics conserves the quantum Hamiltonian.

This is sometimes called the ‘asymmetric-split’ correlation function, since the Boltzmann operator is placed asymmetrically on one side of $\hat{A}$. To picture this function as in Fig. 1a we insert identities into Eq. (18), which when $\hat{A}$ and $\hat{B}$ are functions of position only gives [24]

[TABLE]

We can therefore imagine starting from point $x$ in Fig. 1a and taking an imaginary time path $e^{-\beta\hat{H}}$ ending at $y$, at which $A(y)$ is evaluated. We then take a backwards real time path $e^{i\hat{H}t/\hbar}$ from $y$ to $z$, at which $B(z)$ is evaluated, followed by a real time path $e^{-i\hat{H}t/\hbar}$ from $z$ to $x$, completing the trace.

However, the correlation function is not necessarily real, even for an autocorrelation function (where $\hat{A}=\hat{B}$); one can show by exploiting $[e^{-\beta\hat{H}},e^{\pm i\hat{H}t/\hbar}]=0$ that for arbitrary $\hat{A}$ and $\hat{B}$

[TABLE]

2.2.2 Symmetric-split time-correlation function

Since Eq. (18) can be complex and the classical correlation function is not, we wish to rewrite Eq. (18) to be real. A simple way to do this would be to take the real part of Eq. (18), giving

[TABLE]

which is pictured in Fig. 1b. Although this looks more complex that Eq. (18), if we insert identities as in Eq. (19) and then change to sum-and-difference variables $q=(x+y)/2$, $\Delta=y-x$, noting that the Jacobian of the transformation is unity, we obtain (for $\hat{A}$ which is a linear function of $\hat{x}$)

[TABLE]

We can, for linear operators, consider $\hat{A}$ to be acting at the mid-point of the imaginary time trajectory (this can also hold for some nonlinear operators, see Section 6).

2.2.3 Kubo-transformed time-correlation function

Although Eq. (21) is real and therefore an improvment upon Eq. (18) for classical approximation, the action of $\hat{A}$ at specific points in imaginary time (rather than smoothed over all points) leads to difficulties with classical approximations, as we shall see later. A correlation function which treats all points in imaginary time equally is the Kubo-transformed correlation function [25]

[TABLE]

which corresponds to the zero-time operator $\hat{A}$ being ‘smeared’ through the imaginary time operator $e^{-\beta\hat{H}}$, as pictured in Fig. 1c. This can be obtained for some quantum mechanical properties using linear response theory [26]. In addition to the symmetry properties for Eq. (18), by switching integration limits one can show that the Kubo transformed correlation function is always real,

[TABLE]

and that it obeys detailed balance, i.e.

[TABLE]

and so is more ‘classical’ than the correlation function in Eq. (18). Further symmetry properties of these correlation functions are given in Ref. [27].

2.2.4 Generalized Kubo-transformed time-correlation function

It is possible to rewrite the Kubo-transformed correlation function in a more symmetric form, known as the Generalized Kubo Transformed correlation function [28, 29, 30, 31]. To sketch how this comes about, consider dividing up the imaginary time trajectory $e^{-\beta\hat{H}}$ in the symmetric-split Eq. (22) into $N$ chunks, and at each chunk inserting $e^{-i\hat{H}t/\hbar}e^{i\hat{H}t/\hbar}$, as pictured in Fig. 1d for $N=3$. This gives

[TABLE]

where (for linear $\hat{A}$ and $\hat{B}$)

[TABLE]

with $\hat{A}_{k}$ acting on the $k$th path-integral ‘bead’ and likewise for $\hat{B}$, where we loosely define $q_{i}$ to be the $i$th bead (see appendix A for a discussion of ring polymers and bead terminology). One can show (by evaluating the summations in the correlation function term-by-term and removing $e^{i\hat{H}t/\hbar}e^{-i\hat{H}t/\hbar}=\hat{1}$ identities) that with $\hat{A}$ and $\hat{B}$ defined as in Eq. (27) then this is equal to the conventional Kubo transformed correlation function in the large $N$ limit [32]

[TABLE]

Nonlinear operators [which cannot easily be written as a sum like Eq. (27)] are required for Quantum Transition-State Theory, and are detailed in Section 6. As we shall see later, the advantage of rewriting Eq. (23) as the Generalized Kubo form is that the latter is symmetric with respect to permutation in imaginary time $\tau=\beta_{N}\hbar$, corresponding to permuting the co-ordinates $q_{i}\to q_{i+1}$ [32].

The above is not an exhaustive list of quantum time-correlation functions; there are theoretically infinitely may ways to split the zero-time operator within the Boltzmann distribution[33, 30], one other common technique being $e^{-\beta\hat{H}/2}\hat{A}e^{-\beta\hat{H}/2}$ [33, 24].

By inserting energy eigenstates into Eq. (18) and Eq. (23) one can relate the spectrum of the conventional and Kubo-transformed correlation functions[27, 34]

[TABLE]

where the spectrum is given by

[TABLE]

and likewise for $I_{AB}(\omega)$.

2.2.5 Applications

To illustrate the scope of correlation functions, we now sketch how they may be used to compute diffusion, rates and spectra.

The diffusion constant is obtained as the integral of the Kubo-transformed velocity-velocity autocorrelation function[35]

[TABLE]

where $Z$ is the partition function of the system. The rate constant can be obtained from the long-time limit of the flux-side time-correlation function[36, 33, 26, 37] (of the asymmetric, symmetric Kubo-transformed, and many other forms [33])

[TABLE]

where $Q_{\rm r}(\beta)$ is the partition function in the reactant region and the flux-side correlation function is

[TABLE]

although Eq. (32) also holds for the Kubo-transformed correlation function amongst others [33]. The flux operator is $\hat{F}=[\delta(\hat{q}-q^{\ddagger})\hat{p}+\hat{p}\delta(\hat{q}-q^{\ddagger})]/2m$ where $\delta(x)$ is the Dirac delta function and $q^{\ddagger}$ is the location of the position-space dividing surface. Using the quantum mechanical continuity equation one can show that the exact quantum rate is independent of the location of the dividing surface [38]. The heaviside function $h(\hat{q}-q^{\ddagger})$ is defined such that

[TABLE]

Since the flux operator is the time-derivative of the heaviside operator, the flux-side function is the integral of the flux-flux function [33]

[TABLE]

where $c_{\rm ff}(t)$ is obtained by changing $h(\hat{q}-q^{\ddagger})$ for $\hat{F}$ in Eq. (33), and $c_{\rm fs}(t)$ is minus the derivative of the side-side function

[TABLE]

where $c_{\rm ss}(t)$ is obtained by changing $\hat{F}$ for $h(\hat{q}-q^{\ddagger})$ in Eq. (33). These identities, which generally hold for most classical flux-side time correlation functions too, will prove useful later.

For infra-red spectra, the absorption coefficient is given as[34]

[TABLE]

where $\tilde{I}_{\bm{\mu\mu}}(\omega)$ is the Kubo-transformed dipole autocorrelation function found using Eq. (30), $\mathcal{V}$ corresponds to the volume, $c$ the speed of light and $n(\omega)$ the refraction coefficient (approximately unity in the gas phase).

The above is not exhaustive; other observables can be obtained from thermal quantum time-correlation functions such as neutron scattering [39].

2.3 Moyal series

Having given the exact quantum time-correlation functions in the conventional operator representation, we now consider how the Wigner transform and Moyal series which can be used to rewrite correlation function in terms of phase-space positions and momenta. We use the conventional Kubo-transformed function in this section, but the derivation is equally applicable to the asymmetric or symmetric-split forms.

Inserting position-space identities followed by changing to sum and difference variables as in Eq. (22) gives

[TABLE]

where we have abbreviated the Kubo transform as

[TABLE]

and $\hat{B}(t)=e^{i\hat{H}t/\hbar}\hat{B}e^{-i\hat{H}t/\hbar}$ is the Heisenberg time-evolved $\hat{B}$. We can now insert another identity

[TABLE]

where we have written the Dirac delta function on the first line as its Fourier transform on the second, and convert the $\Delta$ to $-\Delta^{\prime}$ in the second bra-ket of Eq. (40), giving

[TABLE]

where $[\hat{O}]_{\mathrm{W}}$ defines the Wigner transform of operator $\hat{O}$ [40]

[TABLE]

All we have done in Eq. (40)–Eq. (44) is to rewrite the correlation function is terms of classical-like phase-space variables $p$ and $q$. No approximation has been made, an in general solving Eq. (43) exactly is just as difficult as solving the original Eq. (23). The advantage of writing in a classical-like form is the ability to make approximations to the correlation functions such that they can be evaluated using classical or classical-like dynamics.

We now obtain the Liouvillian for a Wigner-transformed correlation function, starting by differentiating Eq. (43) w.r.t. $t$,

[TABLE]

where the commutator arises from noticing $\frac{d}{dt}e^{i\hat{H}t/\hbar}\hat{B}e^{-i\hat{H}t/\hbar}=(i/\hbar)[\hat{H},e^{i\hat{H}t/\hbar}\hat{B}e^{-i\hat{H}t/\hbar}]$. The evaluation of the Wigner transform of the commutator is detailed in Ref. [41] and here we give the main steps.

Using Eq. (17) we can write (dropping the prime on $\Delta^{\prime}$ for simplicity)

[TABLE]

Using the definition $\hat{p}=-i\hbar\frac{d}{d\hat{q}}$, we can take the position derivatives outisde the bra-kets, and using partial differentation show

[TABLE]

and using integration by parts $\frac{d}{d\Delta}$ can be converted into $ip/\hbar$. Combining the above into Eq. (46a) gives

[TABLE]

which is Newton’s first law. For the potential term in Eq. (46b), we observe

[TABLE]

Combining this with $\Delta$ being equivalent to $-i\hbar\frac{d}{dp}$ acting on the entire Wigner Transform we obtain

[TABLE]

where the arrows indicate in which direction the derivative acts, and which is like Newton’s second law with higher-order terms in $\hbar$, as can be seen from expanding the sine series. Combining Eq. (48) and Eq. (50) we obtain

[TABLE]

where $\mathcal{L}_{\rm Moy}$ is the Moyal series[42, 43, 41]

[TABLE]

which is referred to as a series since expanding the sine term gives a series in powers of $\hbar^{2}$. The correlation function is therefore

[TABLE]

In general, computing the action of the Moyal series upon an obserable is as difficult as solving the Schrödinger equation by conventional matrix-based methods, due to the presence of the higher-order derivatives in Eq. (52), although there have been some approaches to address this [44]. In the following sections we therefore explore approximating the Moyal series or generalization of it to obtain classical-like dynamics.

3 LSC-IVR

Arguably the simplest way to approximate $\mathcal{L}_{\rm Moy}$ is to truncate in powers of $\hbar$, giving

[TABLE]

which corresponds to purely classical evolution of the phase-space density from an initial quantum Boltzmann distribution, and has the appealing feature that the error from exact quantum evolution $\mathcal{L}_{\rm Q}$ is known,

[TABLE]

which (by construction) only contains terms of $\mathcal{O}(\hbar^{2})$ and higher. Inserting Eq. (54) into the correlation function gives

[TABLE]

where we have noted that, since $\mathcal{L}_{0}$ is classical, it corresponds to inserting the time-evolved positions and momenta into $[\hat{B}(0)]_{\mathrm{W}}$. Although the Liouvillian has been truncated in powers of $\hbar$, in general this does not mean that the time-evolved observable has been truncated in $\hbar$, since the action of $\frac{\partial}{\partial p}$ in the higher-order terms of $\mathcal{L}_{\rm Moy}$ upon the Wigner transformed obervable ‘brings down’ powers of $\hbar^{-1}$ [45].

The correlation function in Eq. (56) is known as the linearized semiclassical initial value representation (LSC-IVR) or the classical Wigner model, since it can be derived be linearizing the difference in the action between forward-backward trajectories in the semiclassical initial value representation[46], and was later shown to be derivable from linearizing the action of the exact quantum path-integral[47]. The method is exact in the high-temperature limit, for harmonic systems (where the higher terms in the Moyal series vanish without approximation) and as $t\to 0$[32, 47, 46]. LSC-IVR gives fairly good short-time dynamics, though can miss interference effects in non-dissipative systems[48, 6]. A more serious shortcoming is that the classical dynamics does not conserve the quantum Boltzmann distribution, leading to zero-point energy flowing from high-frequency modes to translations and giving spurious effects in simulations[49]; an effect sometimes called ‘zero-point energy leakage’. Evaluating the Wigner-transformed Boltzmann distribution requires a multidimensional Fourier transform which is often approximated [6], and at low temperatures this distribution can have negative values [50]. Nevertheless, it has successfully been applied to reaction rates [51], vibrational energy relaxation and spectra[6, 49].

4 Matsubara dynamics

We have seen how to derive the exact quantum Liouvillian, the Moyal series, and how its truncation to $\mathcal{O}(\hbar^{0})$ leads to classical trajectories, though does not conserve the distribution. This motivates considering whether there are other truncations which give classical trajectories (single derivatives in the Liouvillian) but which also conserve the quantum Boltzmann distribution. Here we show that by truncating in the higher path-integral normal modes a classical, Boltzmann preserving ‘Matsubara’ dynamics is produced. Unfortunately it suffers from the sign problem so is not at present a practical method, though we shall subsequently show how its further approximation leads to the successful approximate methods of CMD, RPMD and TRPMD.

The full derivation of Matsubara Dynamics is in Ref. [32]; here we outline the necessary steps for a one-dimensional system where $\hat{A}$ and $\hat{B}$ are only functions of ${\bf q}$; generalization to more general operators being straightforward[32]. We also require $\hat{A}$ and $\hat{B}$ to be invariant w.r.t. cyclic permutation of the beads $\{q_{i}\}$, which is immediately satisfied if $\hat{A}$ and $\hat{B}$ are linear as in Eq. (27), and is also the case for more general nonlinear operators such as the dividing surface in rate theory[28]. In order to use symmetry w.r.t. imaginary time translation, we use the Generalized Kubo Form in Eq. (26), insert identities and construct a multidimensional Wigner transform as in Eq. (43), giving[32]

[TABLE]

where $N$ is the number of path-integral beads. The Wigner-transformed Boltzmann distribution is given by

[TABLE]

where the bar on $[e^{-\beta\hat{H}}\hat{A}]_{\bar{N}}$ denotes that the bra-kets link together adjacent [$(i-1)$th and $i$th] beads and the real-time evolution is

[TABLE]

where the bra-kets only concern a single bead. As all we have done is insert identities, one could equivalently construct Eq. (57) to have $[e^{-\beta\hat{H}}\hat{A}]_{N}({\bf p},{\bf q})$ and $[\hat{B}(t)]_{\bar{N}}({\bf p},{\bf q})$. However, since the time-evolution bra-kets only concern a single bead, the Liouvillian for Eq. (57) is simply the sum of the Liouvillian in Eq. (52) acting on each bead:

[TABLE]

where

[TABLE]

Truncating Eq. (61) to $\mathcal{O}(\hbar^{2})$ gives LSC-IVR in the same way as truncating $\mathcal{L}_{\rm Moy}$ in Eq. (54) [32].

Formally, one can write the exact correlation function in Eq. (57) as

[TABLE]

although this will generally be even harder to solve exactly than Eq. (43). The benefit of ‘repackaging’ the correlation function as in Eq. (62) is to exploit its symmetry properties w.r.t. imaginary time. For example, $[e^{-\beta\hat{H}}\hat{A}]_{\bar{N}}({\bf p},{\bf q})$ and $[\hat{B}(t)]_{N}({\bf p},{\bf q})$ (as well as the Liouvillian in Eq. (61)) are invariant to cyclic permutation in imaginary time (changing $q_{i}\to q_{i+1}$), whereas this is not obvious with the conventional Kubo-transformed correlation function in Eq. (43). As we shall see later, invariance to translation in imaginary time has a close relationship to the dynamics conserving the quantum Boltzmann distribution.

Instead of writing the correlation function in terms of individual beads, we now consider writing in terms of path-integral normal modes, transforming $({\bf q},{\bf p},\bm{\Delta})\to({\bf Q},{\bf P},{\bf D})$ where the normal modes are numbered $-(N-1)/2\leq j\leq(N-1)/2$ as detailed in Appendix B222Here we consider $N$ and $M$ to be odd for algebraic convenience, even $N$ and $M$ leads to the same result [32].. In brief, the normal modes conventionally originate from diagonalizing the ring-polymer Hamiltonian (see Eq. (159) and Ref. [52]) but here help in evaluating the complex quantum Boltzmann distribution in $[e^{-\beta\hat{H}}\hat{A}]_{\bar{N}}({\bf p},{\bf q})$ and allow an intuitive understanding of the path integral. The lowest mode $Q_{0}$ is (in this definition) the centroid [53, 54, 55], the average position of the beads, and $P_{0}$ the associated momentum. Qualitatively, the modes $Q_{\pm 1}$ describe the size or stretch of the ring polymer [56], $Q_{\pm 2}$ its curvature and so on. $Q_{0}$ can therefore be considered the most ‘classical’ of the modes and the modes are more ‘quantum’ with increasing $|j|$.

In normal modes the correlation function becomes [32]

[TABLE]

where the Liouvillian in normal modes is

[TABLE]

and the potential in normal modes is given by

[TABLE]

If we were to truncate to $\mathcal{O}(\hbar^{0})$ we would recover LSC-IVR once again [32]. Instead, we make a different approximation, truncating from all $N$ to the lowest $M$ path-integral normal modes. From an intuitive perspective, at zero time the highest $N-M$ modes cannot contribute to the (static) correlation function as they are constrained to zero by the quantum Boltzmann operator. One would expect them only to affect the dynamics at longer times when they couple due to anharmonicity in the potential (in a perfectly harmonic potential, the dynamics is separable and the ring polymer normal modes move independently). In the $N\to\infty$ limit, this truncation gives

[TABLE]

and we can therefore define an error Liouvillian [57]

[TABLE]

which is given in full in appendix C.

How many of the lowest $M$ modes should be included? For any physical, analytic potential (one which is smooth, continuous and continuously differentiable) there will be a maximum frequency (second derivative), and provided the frequency of the highest Matsubara mode (see below) is greater than this, all statistical information will be correctly captured (as modes $j\gg M/2$ will move adiabatically to the potential).

For any $M$, the limit $N/M\to\infty$ is taken, and all higher derivatives in Eq. (68) vanish without approximation, since the $l$th derivative scales as $(M/N)^{l-1}$. Consequently333Strictly speaking, the potential in Eq. (68) is $U^{[N]}({\bf Q})$ and this becomes $U^{[M]}({\bf Q})$ after integrating out the non-Matsubara modes detailed below.

[TABLE]

and the single derivatives mean that the dynamics is classical, with a smoothed “Matsubara potential” $U^{[M]}({\bf Q})$.[32]

Because the higher normal modes are not present in the dynamics, nor in $B({\bf Q})$, the higher path-integral momenta can be integrated out from the distribution444This also assumes that $\hat{B}$ is not a function of the higher normal modes in momenta.. This allows the higher-frequency ‘stretch’ variables $\{D_{j},|j|>(M-1)/2\}$ to be integrated out from the distribution. In the $N\to\infty$ limit the Boltzmann bra-kets can be evaluated analytically, leading to the remaining $M$ ${\bf D}$ variables being integrated out by steepest descent. Finally, the higher normal modes in ${\bf Q}$ (which are not affected by $\mathcal{L}^{[M]}$) can be removed by steepest descent. This leads to the classical-like Matsubara correlation function[32]

[TABLE]

where the Matsubara Hamiltonian is

[TABLE]

The phase factor is given by

[TABLE]

where

[TABLE]

are the Matsubara frequencies[58], after which the dynamics is named[32]. Note that, in this definition, the frequencies can be negative since $\tilde{\omega}_{-j}=-\tilde{\omega}_{j}$. $\alpha=\hbar^{1-M}[(M-1)/2]!^{2}$, and the integrals are now implicitly $M$-dimensional as the $N-M$ non-Matsubara modes have been integrated out.

The truncation in normal modes is illustrated pictorially in Fig. 2 and mathematically in Fig. 3.

Since the dynamics in $\mathcal{L}^{[M]}$ is equal to that generated by $H_{M}({\bf P},{\bf Q})$, i.e. $\mathcal{L}^{[M]}=\{\cdot,H_{M}\}$ where $\{\cdot,\cdot\}$ is the Poisson bracket, the dynamics will conserve $H_{M}({\bf P},{\bf Q})$. To show conservation of the phase factor one can either evaluate $\mathcal{L}^{[M]}\theta_{M}({\bf P},{\bf Q})$ and show by trigonometric identities that this vanishes, or use Noether’s theorem[32]. Using the latter method here, we note that the Hamiltonian and therefore the Lagrangian

[TABLE]

is invariant w.r.t. translation in imaginary time. Using straightforward differentiation and that $\frac{d}{dt}Q_{j}=P_{j}/m$,

[TABLE]

and by expanding $\frac{dQ_{j}}{d\tau}$ in bead co-ordinates and applying trigonometric identities we find

[TABLE]

meaning that

[TABLE]

and therefore $\mathcal{L}^{[M]}e^{-\beta[H_{M}({\bf P},{\bf Q})-i\theta_{M}({\bf P},{\bf Q})]}=0$, such that the Matsubara distribution is conserved by the Matsubara Liouvillian, and $C_{AB}^{[M]}(t)$ obeys detailed balance.

Matsubara dynamics is therefore classical and conserves the distribution, but the phase factor in the distribution means that the correlation function is not amenable to computation in large systems. However, for the model systems for which it has been computed, it is more accurate than LSC-IVR, CMD or RPMD[32, 57], and is exact for the position-squared correlation function in a harmonic potential[59] which is not the case for RPMD or CMD[60].

5 Approximations to Matsubara Dynamics

The accuracy of Matsubara dynamics and its intractable nature in large systems suggests that approximations to it which avoid the sign problem may prove more useful in practical applications. Obviously these approximate methods will not in general be as accurate as Matsubara dynamics and one must therefore choose the approximation carefully, in order to remove the sign problem but also keep the dynamics real and preserve the quantum Boltzmann distribution.

In this article we explore three approximations to Matsubara dynamics which fulfil these criteria; a mean-field approximation which yields centroid molecular dynamics (CMD), and moving the momentum contour in the complex distribution of Eq. (69), followed by approximating the resulting complex dynamics deterministically, giving RPMD, or stochastically, giving TRPMD. The full mathematics is given in a series of recent articles [57, 59] and for simplicity only the main details are given here.

5.1 Contour integration

For $t=0$, one can perform contour integration in the complex distribution in Eq. (69), defining

[TABLE]

for all the normal modes. There is no phase factor associated with the centroid ($\tilde{\omega}_{0}=0$), and so the countour of the centroid remains unchanged, which will become important later. Using this transformation, for which the Jacobian is unity, we obtain

[TABLE]

where $R_{M}(\tilde{\bf P},{\bf Q})$ is the ring polymer Hamiltonian in Matsubara modes [57],

[TABLE]

In itself, Eq. (78) is an exact rewriting of Eq. (69), where $\tilde{\bf P}$ are presently complex. However, at zero time, we can evaluate $\{\tilde{P}_{j}\}$ integrals along the real axis, noting that the edges of the contour vanish, giving

[TABLE]

The contour integral is illustrated pictorially in Fig. 4.

At finite time, moving the contour in $\{\tilde{P}_{j}\}$ leads to $\mathcal{L}^{[M]}$ generating complex trajectories which are inherently unstable [61, 62, 63, 64], i.e. we will have exchanged a complex distribution and real dynamics for a real distribution and complex dynamics, and the problem will be equally (if not more) intractable. However, we will see below that moving the contour and discarding (or replacing) undesirable parts of $\mathcal{L}^{[M]}$ can lead to tractable dynamics.

5.2 CMD

If the observables $A({\bf Q})$ and $B({\bf Q})$ are only functions of the centroid $Q_{0}$, we formally rewrite Eq. (69) as

[TABLE]

where the primes denote integration over all modes except $P_{0}$ and $Q_{0}$. We can then define the reduced centroid density

[TABLE]

and differentiation, followed by integration by parts gives

[TABLE]

where the centroid motion alone is given by

[TABLE]

and we have noted that $(\mathcal{L}^{[M]}-\mathcal{L}_{0})\theta_{M}({\bf P},{\bf Q})=0$. At present no approximation has been made and in general direct evaluation of Eq. (83) would be just as difficult as Eq. (69) as the force on the centroid in Eq. (84) requires evaluting the dynamics of all the other normal modes. However, we can define a mean-field force by averaging over all the non-centroid normal modes,

[TABLE]

and then perform contour integration as in Eq. (78) to obtain

[TABLE]

where the normalization is

[TABLE]

We can then approximate the force on the centroid as

[TABLE]

where $F_{\rm f}(Q_{0})$ is defined by Eq. (88), and by discarding $F_{\rm f}(Q_{0})$ we obtain

[TABLE]

from which we can define a centroid-only Liouvillian

[TABLE]

and a formal solution

[TABLE]

We can now perform the contour integration inside $b(Q_{0},P_{0},0)$ giving $b(Q_{0},P_{0},0)=Z_{0}B(Q_{0})$ where $Z_{0}$ is the centroid-density distribution given in Eq. (87). Since $\mathcal{L}_{\rm C}Z_{0}=0$, we can ‘leave’ the distribution at zero time and only propagate $B(Q_{0})$, giving an approximate correlation function

[TABLE]

which is CMD[3, 57, 47, 65, 66, 67, 68, 69]. Consequently, CMD can be obtained from exact quantum dynamics by discarding the motion of the high-frequency modes to obtain Matsubara dynamics, and then making the mean-field approximation $\frac{\partial U^{[M]}({\bf Q})}{\partial Q_{0}}\simeq F_{0}(Q_{0})$, i.e. that the fluctuations around the centroid are negligible. In some situations such as high temperatures this is a reasonable approximation, but at low temperatures where the ring polymer is highly delocalised this can lead to the curvature problem [34] where spectra are artificially broadened and red-shifted, and reaction rates for asymmetric systems are overestimated since the higher normal modes form part of the optimal dividing surface [56]. Because the higher normal modes are integrated out in CMD, it is inaccurate even at $t=0$ for nonlinear operators [60, 70], though various techniques to address this have been proposed [71, 60].

Because $\mathcal{L}_{\rm C}Z_{0}=0$, CMD conserves the distribution function and obeys detailed balance.

In theory, there is no mathematical obligation to take the mean field of all non-centroid modes, and one could average out over a subset, such as the most highly oscillatory ones. While this would include some level of fluctuations, the distribution of the non-centroid modes which were not integrated out would still suffer from the sign problem.

5.3 RPMD

As noted in section 5.1, analytic continuation of the non-centroid momenta is mathematically possible, and the integrand can be proven to be holomorphic in that region of the complex plane[59], meaning that there are no singularities to worry about. The complex Liouvillian can be written as its real and imaginary parts,

[TABLE]

where

[TABLE]

is the ring polymer Liouvillian (using Matsubara frequencies) and

[TABLE]

One can show that both $\mathcal{L}^{[M]}_{\Re}$ and $i\mathcal{L}^{[M]}_{\Im}$ separately conserve the distribution in Eq. (80), and so discarding $i\mathcal{L}^{[M]}_{\Im}$ leads to a correlation function with a real distribution and a real dynamics which conserves it,

[TABLE]

which is RPMD[57, 27]. This means that the error in the evolution between exact quantum dynamics and RPMD can be stated in closed form as the error between exact quantum dynamics and Matsubara dynamics [Eq. (160)], followed by a contour integral and discarding $\mathcal{L}^{[M]}_{\Im}$ [Eq. (96)]555Strictly speaking, one also discards the vertical edges of the integral contour, which are believed to be zero[59]..

Since $\mathcal{L}^{[M]}_{\rm RP}e^{-\beta R_{M}(\bar{\bf P},{\bf Q})}=0$, RPMD conserves the distribution and $C_{AB}^{\rm RP}(t)$ obeys detailed balance. Strictly speaking, Eq. (97) is RPMD with Matsubara frequencies, but in the $M\to\infty$ and $N/M\to\infty$ limits (implicitly taken here), only the lowest Matsubara modes will participate in the statistics and dynamics, the others being constrained to zero by the spring terms in $R_{M}(\bar{\bf P},{\bf Q})$, and correlation functions employing Matsubara and ring polymer frequencies will converge to the same result [57].

One unfortunate effect of discarding $\mathcal{L}^{[M]}_{\Im}$ is that it shifts the frequencies of the non-centroid normal modes; in a harmonic potential $V(q)=\frac{1}{2}m\omega_{h}^{2}q^{2}$, they become [4]

[TABLE]

This leads to the so-called ‘spurious resonances’ problem in spectra, where resonances between ring polymer frequencies and physical frequencies (such as stretching vibrations) lead to spurious extra spectra peaks which are temperature-dependent[70, 34, 72, 73].

5.4 TRPMD

To address the artificial shifting of frequencies upon discarding $i\mathcal{L}^{[M]}_{\Im}$, we consider replacing it with an operator which will conserve the distribution but also provide the correct oscillation frequency. The standard analysis of a damped harmonic oscillator[2] shows that a friction term will reduce the oscillation frequency, so we consider defining[59]

[TABLE]

where $\mathcal{A}_{\rm wn}^{[M]{\dagger}}$ is the adjoint of a white-noise Fokker-Planck operator[1],

[TABLE]

The first term on the RHS of Eq. (100) corresponds to the drag cause by the semidefinite friction matrix $\bm{\Gamma}$ (which we assume is diagonal in what follows) and the second term represents the ‘kicks’ imparted to the individual momenta of stochastic trajectories[2]. Inserting Eq. (99) into the analytically continued correlation function gives

[TABLE]

which is TRPMD[59, 70]. Similar to RPMD, the approximation in the dynamics between exact quantum evolution and TRPMD is therefore known, namely $\mathcal{L}_{\rm er}$ followed by a contour integral and replacing $i\mathcal{L}^{[M]}_{\Im}$ with $\mathcal{A}_{\rm wn}^{[M]{\dagger}}$.

Using integrating by parts one can obtain the (non-adjoint) of the Fokker-Planck operator in Eq. (100) as [1]

[TABLE]

such that the Eq. (101) can be rewritten as

[TABLE]

We can then show that $\mathcal{A}_{\rm RP}^{[M]}e^{-\beta R_{M}(\bar{\bf P},{\bf Q})}=0$ such that the stochastic dynamics of the system conserves the distribution. Showing that the correlation function obeys detailed balance is more complicated (since $\mathcal{A}_{\rm RP}^{[M]}$ contains double derivatives) and this is detailed in Ref. [23].

Defining the friction matrix to be $\mathbf{\Gamma}_{jk}=2|\tilde{\omega}_{j}|\delta_{jk}$ leads to the correct oscillation frequency of all ring polymer normal modes in a harmonic potential, and therefore give the correct zero-time value and oscillation frequency for the harmonic position-squared autocorrelation function [59], which neither RPMD nor CMD can achieve [60, 59]. More importantly for spectra, a friction matrix of $\mathbf{\Gamma}_{jk}=\sqrt{2}|\tilde{\omega}_{j}|\delta_{jk}$ will lead all peaks in the position autocorrelation function for a harmonic oscillator to be at the correct (external) frequency, and therefore provides a unique value of $\mathbf{\Gamma}_{jk}$ for computation of spectra which is between the values previously suggested on the basis of optimal sampling[52, 70].

Although TRPMD improves on both CMD and RPMD for spectra [70], the friction causes unphysical slowing of reaction rates beneath the crossover temperature [23].

5.5 Summary

The various approximations used to obtain LSC-IVR, CMD, RPMD and TRPMD are illustrated schematically in Fig. 5 and their properties summarized in Table 1. For many systems with mild quantum effects some or all of these methods will produce similar results [74], and all are exact in the high-temperature (classical) limit [70, 27, 6, 68], the $t\to 0$ limit [75, 70, 6] and for the position autocorrelation function of a harmonic oscillator [6, 70, 68, 27, 57]. Although we have shown that CMD can be obtained directly from Matsubara dynamics as a mean field approximation, it can also be obtained as a mean field approximation to RPMD and TRPMD using the same methodology, as shown for RPMD in Ref. [76].

6 Quantum transition-state theory

Having considered time-correlation functions, we now consider one of their principal applications: reaction rate calculation, and how the foregoing mathematical ‘toolkit’ can be used to obtain quantum transition-state theory.

6.1 Background

Here we provide a brief outline of the development of rate theory to place the material discussed here in context; for a fuller historical overview see Ref. [80].

The earliest widely-accepted rate formula is arguably the Arrhenius equation

[TABLE]

where $A$ is the pre-exponential (frequency) factor and $E_{a}$ is the activation energy. Obtained empirically, there was originally no clear prescription for determining $A$ a priori. In 1935 Eyring[81, 82] along with Evans and Polanyi[83] proposed

[TABLE]

where $\sqrt{\frac{m}{2\pi\beta\hbar^{2}}}K^{*}$ is the equilibrium constant between the reactants and the activated complex (the thermal probability of finding the system at the transition state), $1/\sqrt{2\pi m\beta}$ is the thermal flux and, to quote Eyring[82]

The transmission coefficient $\kappa$ is just the ratio of systems crossing the barrier to systems reacting…Fortunately, as stated for many reactions we make a negligible error by taking it as unity.

Consequently, Eq. (106) (hereafter “Eyring TST”) is the thermal flux multiplied by the probability of forming the activated complex, or in modern terminology, the thermal flux through the dividing surface, which gives the exact rate if there is no recrossing. The partition functions involved are calculated quantum mechanically, but the motion through the transition state is assumed to be classical and separable from motion orthogonal to the dividing surface, which is not always the case [84] and in some circumstances can lead to considerable errors.

6.2 Classical rate theory

Determining the functional form of the transmission coefficient was placed on a firmer theoretical footing in the 1970s by constructing a classical flux-side correlation function to determine the classical rate[85, 86],

[TABLE]

where (in one dimension for simplicity)

[TABLE]

This correlates the flux through $q^{\ddagger}$ at zero time, $\delta(q-q^{{\ddagger}})p/m$, with whether the system is in the product region at time $t$, $h(q_{t}-q^{\ddagger})$. Here $Q_{\rm r}(\beta)$ is the partition function in the reactant region, $\delta(q-q^{\ddagger})$ is a Dirac delta function and $h(q_{t}-q^{\ddagger})$ a heaviside function, similar to the quantum case. For an $F$-dimensional system one defines a reaction co-ordinate $f({\bf q})$ such that $f({\bf q})=0$ defines an $(F-1)$-dimensional dividing surface, $f({\bf q})>0$ is the product region and $f({\bf q})<0$ is the reactant region.

Strictly speaking, the infinite-time limit in Eq. (107) is only valid for gas-phase scattering. For condensed-phase systems, in order to define a rate there must be sufficient separation in timescales between reaction and equilibration for plateau in $c_{\rm fs}(t)$ to emerge, at which point the rate is evaluated[85].

6.3 Classical TST

Here we show how the classical TST rate is related to the short-time limit of Eq. (108) and therefore to the classical rate. In the process we obtain an algebraic expression for the transmission coefficient. We firstly formally rewrite Eq. (108) as

[TABLE]

where $\mathcal{L}$ is the classical Liouvillian given in Eq. (7), and we have used the algebra in Section 2 to take $e^{\mathcal{L}t}$ ‘inside’ the heaviside function, since $\mathcal{L}$ only contains single derivatives in $p$ and $q$. Because the heaviside function is discontinuous, one has to be careful expanding $e^{\mathcal{L}t}h(q-q^{\ddagger})$ around $t=0$, and it is mathematically simpler to use Eq. (109b) rather than Eq. (109a).

In the short-time limit,

[TABLE]

We then note that the Dirac delta function constrains $q=q^{\ddagger}$ and that the heaviside function is invariant to the scaling of its argument, such that

[TABLE]

Putting Eq. (111) back into Eq. (108) gives

[TABLE]

where the integrals in $p$ and $q$ have become separable. The momentum integral is proportional to the thermal flux at inverse temperature $\beta$, and the position integral is proportional to the thermal probability of reaching the transition state $q^{\ddagger}$. Comparing this with Eq. (106), we see that this (suitably scaled by the partition function $Q_{\rm r}(\beta)$) is the classical transition-state theory rate,

[TABLE]

The transmission coefficient, which is the ratio of the classical TST rate to the exact classical rate, is therefore given by

[TABLE]

where

[TABLE]

In practice, rates are often calculated using expressions such as Eq. (115), known as the Bennett-Chandler factorization [21], since this splits the calculation into a statistical part $k_{\rm cl}^{\ddagger}(\beta)$ for which there exists a huge repertoire of efficient sampling techniques[21, 24], and a dynamical part $\kappa(t)$ which can be obtained from a molecular dynamics simulation.

From this we can also obtain a mathematical criterion for recrossing. We firstly note that from Eq. (114), $\lim_{t\to 0_{+}}\kappa(t)=1$, and obtain the time-derivative of $\kappa(t)$ (c.f. Eq. (37)),

[TABLE]

where the classical flux-flux correlation function is

[TABLE]

This gives the flux of particles through the barrier at time $t$, which also went past the barrier at time $t=0$, i.e. the extent of recrossing. If there is no recrossing then $c_{\rm ff}(t)=0$ for all $t>0_{+}$, $\kappa(t)=1$ for all $t\geq 0$, and $k_{\rm cl}(\beta)=k_{\rm cl}^{\ddagger}(\beta)$ which fulfils Eyring’s requirement for a TST.

We can therefore mathematically define classical TST as a rate theory fulfilling two simple criteria:

$k_{\rm cl}^{\ddagger}(\beta)=\frac{1}{Q_{\rm r}(\beta)}\lim_{t\to 0_{+}}c_{\rm fs}(t)$ such that 2. 2.

$k_{\rm cl}^{\ddagger}(\beta)=k_{\rm cl}(\beta)$ if $c_{\rm ff}(t)=0$ for all $t>0_{+}$.

These criteria are not new and are essentially a mathematical summary of the generally-accepted definition of classical transition-state theory [87, 85, 7, 88, 21].

We now briefly note further properties of classical TST which will be useful to compare to QTST. First, if the flux-side time correlation function was defined with two dividing surfaces in different places

[TABLE]

where $q^{\ddagger}_{1}\neq q^{\ddagger}_{2}$ then

[TABLE]

such that

[TABLE]

since the integral in momentum is odd. The existence of a nonzero TST is therefore a consequence of the two dividing surfaces being in the same place [28].

Second, the separability of the position and momentum terms in the classical TST expression Eq. (112b) means that momentum can be integrated out which (along with evaluating the partition function for a scattering system) gives

[TABLE]

showing that classical TST does not require the simultaneous specification of position and momentum, even though this is allowed in classical mechanics.

Third, classical rate theory is independent of the location of the dividing surface [36, 38], which can be shown algebraically by differentiating $c_{\rm fs}(t)$ w.r.t. $q^{\ddagger}$, rearranging, and showing that this corresponds to the system traversing the barrier at time $t$ having starting at the barrier at $t=0$, which cannot be the case at long times if there is a plateau in $c_{\rm fs}(t)$ and the rate is defined. However, classical TST is exponentially sensitive to the dividing surface. Since recrossing only reduces the rate (by the heaviside function discarding trajectories with positive momentum, or including trajectories with initially negative momentum), classical TST is an upper bound to the classical rate. This property can be used to variationally optimize the location of the dividing surface in multidimensional systems [87], since in an $F$-dimensional system the dividing surface is an $(F-1)$-dimensional hypersurface, and locating the position of the optimal dividing surface [the one which minimises $k_{\rm cl}^{\ddagger}(\beta)$ and maximises $\kappa(t)$] is difficult.

In summary, classical transition-state theory is the instantaneous thermal classical flux through a position-space dividing surface, which is equal to the exact (classical) rate in the absence of recrossing ($c_{\rm ff}(t>0)=0$) by the classical dynamics of the system. It also implicitly assumes that the reactants are in thermal equilibrium (and in equilibrium with the transition state) and that the reaction is electronically adiabatic, proceeding on a single Born-Oppenheimer potential energy surface.[7] The advantages of classical TST over full classical rate calculation is computational simplicity, only requiring knowledge of the PES at the dividing surface and no dynamics, and that it is generally easy to tell in advance if TST will provide a good approximation to the rate. TST works for direct reactions where there is a significant thermal barrier between reactants and products (significantly greater than $k_{\rm B}T$); although it is only exact in a small number of cases (such as one dimensional systems with the optimal dividing surface), recrossing of the optimal dividing surface is often small and it is therefore a good approximation, and upper bound, to the rate.[7, 80] It is not expected to work where reactions are diffusive (involving multiple recrossings and therefore a low transmission coefficient), systems with long-lived intermediates (where defining a dividing surface is problematic) or systems with pronounced quantum effects.

6.4 Quantum TST

While very successful for heavy atoms at high temperatures, classical TST does not include any quantum mechanical effects such as tunnelling and zero-point energy, which can lead to significant (many orders of magnitude) deviation between the classical result and the experimental or the quantum result, particularly at low temperatures (see e.g. Ref. [89]). One can, of course, try to include quantum effects into classical TST [7], such as in the standard Wigner-Eyring model where partition functions in modes orthogonal to the reaction co-ordinate are evaluated quantum mechanically, but motion through the saddle point is assumed to be classical and separable to motion orthogonal to it, which is frequently not the case [84].

There is considerable historical debate on the existence of quantum transition-state theory, for which the reader is referred to (for example) Refs. [90, 91, 92, 36, 55, 93, 94, 95]. In short, in the late 1930s Wigner and others considered incorporating quantum effects such as tunnelling into transition-state theory, and noted that there were difficulties due to (a) the non-locality of the quantum Boltzmann operator and (b) the uncertainty principle.

The non-locality of the quantum Boltzmann operator means that the dividing surface must act on a point or points of the imaginary time trajectory embodied in $e^{-\beta\hat{H}}$. The development of path-integral techniques by Feynmann [96] and many others means that the dividing surface can be written as a function of path-integral space, $f({\bf q})$, taking the positions of path-integral beads $q_{1},q_{2},\ldots,q_{N}$ as its argument, such that $f({\bf q})=0$ at the dividing surface. To define a rigorous QTST where the only assumption is no recrossing [93] we therefore have to consider recrossing of the path-integral dividing surface $f({\bf q})$, and recrossing of any surfaces orthogonal to it in path-integral space, which we denote $g({\bf q})$ [29]666Orthogonality formally means than $f({\bf q})\overleftarrow{\nabla}\cdot\overrightarrow{\nabla}g({\bf q})=0$ where $\nabla g({\bf q})$ is the gradient of $g({\bf q})$ [29]..

Concerning the uncertainty principle, specifying the dividing surface in path-integral space allows for a delocalised imaginary-time trajectory and therefore uncertainty in the individual bead positions. We also note that there is no requirement for simultaneous specification of position and momentum in classical TST (see above) and there is no a priori reason why this should be required in the quantum case either.

Extending the definition of classical TST to the quantum case, quantum transition-state theory is therefore defined as the instantaneous thermal flux through a position-dependent dividing surface which gives the exact quantum rate in the absence of recrossing, both of the dividing surface and of the surfaces orthogonal to it in path-integral space[28, 29, 30, 31, 97]. Mathematically, we denote $C_{\rm fs}(t)$ to denote a flux-side function correlating flux through $f({\bf q})$ at $t=0$ with time-evolved side through $f({\bf q})$ (similarly for $C_{\rm ff}(t)$) and $M_{\rm fs}(t)$ for a flux-side function correlating flux through $f({\bf q})$ with time-evolved side through $g({\bf q})$ [29]. The criteria for QTST given algebraically are therefore

$k^{{\ddagger}}_{\rm Q}(\beta)=\lim_{t\to 0_{+}}C_{\rm fs}(t)/Q_{\rm r}(\beta)$ such that 2. 2.

$k_{\rm Q}^{\ddagger}(\beta)=k_{\rm Q}(\beta)$ if $C_{\rm ff}(t)=0$ and $M_{\rm ff}(t)=0$ for all $t>0_{+}$ and all $g({\bf q})$.

We stress that the dynamics in these quantum correlation functions is the exact quantum dynamics ($e^{-i\hat{H}t/\hbar}$) and not any of the approximate quantum methods discussed above.

The historical difficulties of formulating a rigorous QTST (satisfying both of the above criteria) led to the development of a huge range of heuristic quantum mechanical rate theories that used transition-state arguments [53, 54, 55, 98, 36, 99] in addition to alternative approaches such as instanton theory [56, 100, 101], quantum instanton methods [102] and many others discussed elsewhere [9, 8]. There have also been other, generally broader, definitions of QTST in (for example) Refs. [103, 8]. The definition of QTST used in this article is based on Eyring’s original definition of TST and means that one has a priori knowledge of its applicability: provided there is minimal recrossing QTST will be a good approximation to the rate.

6.4.1 Wigner-Miller TST

Having defined QTST we show how to derive a simple expression satisfying the criteria for a QTST, but which is unreliable at low temperatures. In the followed sections we will extend this to obtain an expression which has positive definite Boltzmann statistics, i.e. is guaranteed to be positive at any finite temperature. The original QTST derivation evaluated time-evolution bra-kets algebraically [28]; here we rederive these expressions in the Moyal series formalism, which is arguably simpler.

As in classical mechanics, the key ingredient in formulating a QTST is ensuring that the two dividing surfaces are located in the same place in path-integral space, such that they coalesce in the $t\to 0_{+}$ limit. This has to be done carefully, since the quantum Boltzmann operators is nonlocal, unlike the classical Boltzmann operator. We start with the Wigner-transformed side-side correlation function

[TABLE]

at $t=0$, the dividing surfaces in Eq. (122) are clearly the function of the same co-ordinate and in the same place (they are not separated by an imaginary-time trajectory). We obtain the flux-side correlation function as

[TABLE]

where we have noted that the adjoint of the Liouvillian is its negative[104], and that $\mathcal{L}_{\rm Moy}[e^{-\beta\hat{H}}]_{W}(p,q)=0$ since exact quantum dynamics conserves the quantum Boltzmann distribution. We illustrate Eq. (123) schematically in Fig. 6.

Expanding $e^{\mathcal{L}_{\rm Moy}t}$ in a Taylor series to find the $t\to 0_{+}$ limit is mathematically problematic since $h(q-q^{\ddagger})$ is discontinuous around $q=q^{\ddagger}$, as for the classical case. However, we can instead write

[TABLE]

where $\mathcal{L}_{0}$ is the classical Liouvillian is defined in Eq. (54) and $\mathcal{L}_{\rm Q}$ defined in Eq. (55) contains the higher-order quantum terms. Because $\mathcal{L}_{0}$ only contains single derivatives we can use the maths as for the classical case to show

[TABLE]

and therefore

[TABLE]

Inserting Eq. (126) into Eq. (122) immediately gives

[TABLE]

This is identical to a rate expression introduced heuristically by Wigner in 1932 [99] and was subsequently reintroduced and developed for the description of quantum mechanical reaction rates [105, 50].

The proof that this gives the exact rate in the absence of recrossing is given in [30], fulfilling the second criterion for a QTST. In brief, since the dividing surface acts only on one point in path-integral space (the average of the end-points of the imaginary time path, see Fig. 6), there are no orthogonal surfaces whose recrossing need be considered. Consequently, as the first criterion for QTST is satisfied, one can combine this with Eq. (37) to rewrite the second criterion as $\lim_{t\to\infty}C_{\rm fs}^{[1]}(t)/Q_{\rm r}(\beta)=k_{\rm Q}(\beta)$. This is then proven by evaluating both sides of the equation using quantum scattering theory [30, 36, 106] where the RHS is given by Eq. (32).

While providing a reasonable description at relatively high temperatures, beneath the ‘crossover temperature’ into deep tunnelling (see appendix D) the thermal Wigner distribution becomes non-positive definite, such that Eq. (127) can produce spurious negative rates.[50, 28] This is because only the average of the forward and backward imaginary time paths are constrained to be at the barrier, and the resulting path-integral ‘string’ will sag over the barrier at low temperatures [50, 28].

6.4.2 Positive-definite statistics

To ensure that the rate is positive at any finite temperature, the Generalized Kubo correlation function can be used. The full derivation is given in Refs. [28, 29] and here we sketch the pertinent details. A key part of this is defining a dividing surface in path-integral space $f({\bf q})$ which must separate the products and reactants, converge with $N$ and (in order to maximise the free energy) a permutationally-invariant function of the path-integral beads [28]. In the terminology of Matsubara dynamics, this means that it must be composed of a finite number of $\mathcal{K}$ Matsubara modes.[97]

We start with the Kubo-transformed side-side correlation function

[TABLE]

We then transform the correlation function to path-integral normal modes, without truncating the non-Matsubara modes:

[TABLE]

As before, we differentiate w.r.t. $t$ to obtain the flux-side correlation function

[TABLE]

where $S({\bf P},{\bf Q})$ is the ring-polymer flux

[TABLE]

that is only a function of the lowest $\mathcal{K}$ normal modes. Equation (130) and its short-time limit is given schematically in Fig. 7.

In the short-time limit we can separate the propagator

[TABLE]

where $\mathcal{L}^{[M]}$ is given in Eq. (68) and $\mathcal{L}_{\rm er}=\mathcal{L}_{\rm Moy}^{[N]}-\mathcal{L}^{[M]}$, given in appendix C. For this derivation, we can choose any $M\geq\mathcal{K}$. Using similar algebra to the classical and Wigner-Miller TST cases, we then show

[TABLE]

where we have Taylor-expanded $f({\bf Q})$ and noted that $\mathcal{L}_{\rm er}h[f({\bf Q})+S({\bf P},{\bf Q})t]=0$ since $h[f({\bf Q})+S({\bf P},{\bf Q})t]$ only contains Matsubara modes and all terms in $\mathcal{L}_{\rm er}$ contain derivatives of non-Matsubara modes. This gives

[TABLE]

and inserting Eq. (134) into Eq. (130) we obtain

[TABLE]

which is a nonzero $t\to 0_{+}$ quantum transition-state theory by the first criterion, from which we define $k_{\rm Q}^{\ddagger}(\beta)=\lim_{t\to 0_{+}}C_{\rm fs}^{[N]}(t)/Q_{\rm r}(\beta)$.

To evaluate Eq. (135) we can, without approximation, integrate out the non-Matsubara ${\bf P}$, followed by ${\bf D}$ inside $[e^{-\beta\hat{H}}]_{\bar{N}}({\bf P},{\bf Q})$ and the non-Matsubara ${\bf Q}$ (which by construction are not required to evaluate the distribution) to give

[TABLE]

This expression is identical to the short-time limit of the Matsubara flux-side time-correlation function, or ‘Matsubara transition-state theory’ (M-TST).

To address the phase factor, we then move the contour in ${\bf P}$ to generate a ring polymer potential. If the dividing surface contains non-centroid modes we obtain

[TABLE]

which appears complex, but the imaginary part corresponds to the change in dividing surface with imaginary time $\tau$, which is zero by construction:

[TABLE]

where we have used $\tilde{\omega}_{j}Q_{-j}=-\frac{dQ_{j}}{d\tau}$ from Ref. [32]. This leads immediately to

[TABLE]

which is RPMD-TST with Matsubara frequencies. As for other static and dynamical properties, this is formally identical to RPMD-TST with ring-polymer frequencies in the large $M$, $N\to\infty$ limit considered here. [32]

We have therefore shown that $\lim_{t\to 0_{+}}C_{\rm fs}^{[N]}(t)/Q_{\rm r}(\beta)$ is nonzero giving a QTST by the first criterion. To show that it fulfils the second criterion, we apply Eq. (37) to the second criterion, and note that $\lim_{t\to 0_{+}}M_{\rm fs}(t)=0$ since the dividing surfaces are in different locations in path-integral space. It then becomes sufficient to prove that $\lim_{t\to\infty}C_{\rm fs}^{[N]}(t)/Q_{\rm r}(\beta)=k_{\rm Q}(\beta)$ when $\lim_{t\to\infty}M_{\rm fs}(t)=0$. The mathematics is given in Ref. [29], and in brief the long-time limits are evaluated using quantum scattering theory and we then show that if $\lim_{t\to\infty}M_{\rm fs}(t)=0$ for all $g({\bf q})$ orthogonal to $f({\bf q})$ then $\lim_{t\to\infty}C_{\rm fs}^{[N]}(t)$ is equivalent to the long-time limit of $c_{\rm fs}(t)$ in Eq. (33) which by Eq. (32) fulfils the second criterion.

In theory, it is possible to systematically improve QTST to the exact quantum result by computing the recrossing in $C_{\rm fs}^{[N]}(t)$ and $M_{\rm fs}(t)$ [29, 31], but in practice this is more expensive than a conventional quantum calculation.

6.4.3 Summary

We have rederived RPMD-TST and M-TST from a quantum flux-side time-correlation function using the Liouvillian formalism, finding that both are true quantum transition-state theories. Interestingly, for Matsubara TST to be equivalent to QTST only requires that the dividing surface is a function of a finite number of Matsubara modes, but showing the equivalence to RPMD-TST requires the extra condition that the dividing surface is invariant to cyclic permutation.

We also observe that, when the centroid dividing surface is used, RPMD-TST reduces to the earlier centroid-TST [53, 54, 55, 28]. In fact, a recent article claimed to have derived QTST and found that this was equal to Centroid-TST and not RPMD-TST [107], and which was shown to be an artifact of Ref. [107] only considering a centroid dividing surface [97].

In practice, locating the optimal dividing surface $f({\bf q})$ is complicated and, particularly at low temperatures, may take on a complicated curvilinear form [56]. Because RPMD rate theory is independent of the location of the dividing surface [38], the RPMD rate will be equal to the exact quantum rate is there is no recrossing of the optimal dividing surface (the one which minimises $k_{\rm QM}^{\ddagger}$) or those orthogonal to it in path-integral space by either the exact quantum dynamics or the RPMD dynamics of the system. As for classical TST, in general there will be some recrossing, and consequently RPMD is expected to be a good approximation to the rate.

RPMD rate theory itself has seen a huge range of applications, many of which are discussed in Refs. [4, 5]. To mention a few, after initial application to model systems [88, 38] it was applied to proton transfer [108], bimolecular reaction rates [109, 110] and diffusion in ice and clathrates [111, 112]. QTST has also been applied to improve standard tunnelling corrections [113].

Whereas classical TST is an upper bound to the classical rate, QTST is not a strict upper bound to the quantum rate[28]. However, in general QTST is a good approximation to an upper bound provided that there are not significant coherences in the reaction dynamics [28].

7 Future directions

Having surveyed how CMD, RPMD and TRPMD can be considered as approximations to Matsubara dynamics, we briefly consider areas for further development of the field.

7.1 Nonadiabatic systems

For small or model systems, exact methods can be applied such as MCTDH [114], and the past few decades have seen considerably development of approximate methods. There exist a wide variety of methods to model non-adiabatic processes using classical-like trajectories, including surface-hopping [115, 116, 117], various linearized methods [118], and mixed quantum-classical [119, 120, 121] methods. A common and successful method to map discrete electronic states to continuous classical variables is to use ‘mapping variables’, where singly excited oscillator states are inserted and electronic states represented by their fictitious positions and momenta [122, 123, 124, 125]. There are, of course, many other possible mappings [125] but the simplicity and ease of implementation of mapping variables appears to have led to their widespread application to semiclassical [126, 127], quasiclassical [128], (partially) linearized [129, 130, 131, 18, 132, 133, 134, 135], and path integral dynamics [136, 137, 138]. Although the propagator (Moyal series) for a single surface systems was obtained in 1949 [43], the analogue of this in mapping variables was not derived until 2016 [104].

Despite this progress there remains, to the author’s knowledge, no method which has classical-like scaling in all degrees of freedom, conserves the quantum Boltzmann distribution and reproduces Rabi oscillations, though there are a number of methods which incorporate some of these desirable properties[139]. There is also, at present, no widely-accepted ‘true’ ($t\to 0_{+}$) non-adiabatic quantum transition-state theory with a dividing surface in electronic space—though this does not mean that one does not exist. For a non-adiabatic system with a dividing surface solely in position space, QTST is simply RPMD-TST with a mean-field non-adiabatic potential [31], which means that mean-field non-adiabatic RPMD [140, 141] will provide a good approximation to the exact quantum rate when there is minimal recrossing of the position-space dividing surface by either the (mean field) ring polymer dynamics or the exact quantum dynamics. While this appears to be true for some model systems with large non-adiabatic coupling [140], this is unlikely to hold in regimes of small coupling [141]. Even within existing methods, such as non-adiabatic RPMD, there are a variety of implementations [136, 137, 142, 140, 141] and it is not always clear which one will be superior in any given situation.

7.2 Theoretical development

There may also be the possibility of applying Matsubara dynamics (or a similar approximate quantum dynamics) to the computation of nonlinear response functions [143] which can diverge in a purely classical calculation [144]. There may also be other classical-like approximations to quantum dynamics (and maybe Matsubara dynamics) that for some systems are more accurate[145]. Very recent research has obtained out-of-equilibrium RPMD and CMD from Matsubara dynamics [20], which should be useful tools for excited state quantum dynamics.

7.3 Computational development

For a method to bridge the gap between theoretical development and routine application in large chemical systems, the speed of computation needs to be comparable to that of a standard classical molecular dynamics simulation. There have consequently been a large range of methods developed to implement the approximate methods described here accurately and efficiently.

For single-surface systems, there have been impressive applications including a study of dynamics and dissipation in enzyme catalysis [146] and proton transport in water nanowires [79], though applications to large systems are often limited by the cost of the potential. Various techniques have evolved to address this, including ring polymer contraction [111, 147, 148] and thermostatting [52, 70, 59].

Open source codes such as i-Pi [149] and RPMDrate [150] have been developed to facilitate application to wide-ranging systems.

8 Conclusions

In this New View we have reviewed how a number of successful approximate quantum dynamics methods can be obtained from exact quantum time evolution and used the Liouvillian and Moyal formalisms to rederive quantum transition-state theory.

We have mainly considered the mathematical basis for these theories and shown what terms they discard from the exact quantum evolution to obtain a classical-like dynamics from which to compute a correlation function. Provided the discarded error terms are small, the approximate correlation function will be a good approximation to the exact quantum correlation function. By considering cases where this is (and is not) the case, we can propose a priori situations where a particular methods is likely to work, and therefore advise the usage of approximate methods, summarized in Table 1.

We then revisit classical and quantum transition-state theory and derive QTST in the Matsubara formalism, showing that Matsubara-TST is a true QTST. Provided the dividing surface is permutationally invariant, RPMD-TST is equivalent to Matsubara-TST, unlike the dynamics in RPMD which is only an approximation to Matsubara dynamics. While of limited computational importance by itself (due to the phase factor in the Matsubara distribution) this may facilitate the derivation of other (possibly more accurate) rate theories.

While there has been much progress in recent years, there remain many avenues for further theoretical development. There is arguably no clear consensus on how to apply approximate path-integral methods to non-adiabatic systems, nor a $t\to 0_{+}$ non-adiabatic QTST, the existence of which is an open question. There is also scope for applying the approximate methods discussed here to out-of-equilibrium systems and nonlinear response functions, in addition to developing efficient computational algorithms for implementing these methods in code libraries and for large systems.

9 Acknowledgements

TJHH wishes to thank Jesus College, Cambridge for funding, Stuart Althorpe for helpful discussions, and Srinath Ranya and Elliot C. Eklund for comments on the manuscript.

Appendix A Ring Polymers

There exists a vast literature on ring polymers [96, 151] and here we give the standard derivation of the expression for a partition function [24] for the benefit of those unfamiliar or new to the subject.

For a quantum mechanical partition function

[TABLE]

we can perform the Trotter discretization

[TABLE]

and in the $N\to\infty$ limit, expand $e^{-\beta_{N}\hat{H}}$ symmetrically as

[TABLE]

where $\hat{V}=V(\hat{q})$ and $\hat{T}=\hat{p}^{2}/2m$. We then insert $N$ sets of position identities, $\int dq_{i}|q_{i}\rangle\langle q_{i}|$, $i=1,\ldots,N$,

[TABLE]

where we have noted $e^{-\beta_{N}\hat{V}/2}|q_{i}\rangle=|q_{i}\rangle e^{-\beta_{N}V(q_{i})/2}$ and cyclic permutation within indices to go from Eq. (143a) to Eq. (143b). By inserting momentum eigenstates, we then evaluate

[TABLE]

by contour integration, and by inserting Eq. (144) into Eq. (143b) obtain

[TABLE]

where the ring polymer potential is

[TABLE]

One can re-insert $N$ momentum identities[152]

[TABLE]

in $p_{i}$, $i=1,\ldots,N$ to give

[TABLE]

where the ring polymer Hamiltonian is

[TABLE]

The above derivation is exact for static properties and the dynamics generated by Eq. (149) was originally proposed as a sampling tool[152]. The $\{q_{i}\}$ are known as ring polymer ‘beads’ and in practice their number $N$ is treated as a convergence parameter in a numerical simulation.

Appendix B Normal modes

The ring-polymer normal modes are defined here as in Ref. [57],

[TABLE]

where $j=-N/2+1,\ldots,0,\ldots,N/2$ and likewise for ${\bf P}$, where

[TABLE]

where the $j=N/2$ mode is omitted if $N$ is odd. The transformation is not unitary, but defined such that the normal modes converge in the $N\to\infty$ limit. This leads to frequencies in the complex Boltzmann distribution of

[TABLE]

which, for large $N$ and finite $j$, become the Matsubara frequencies [58]

[TABLE]

The observables $A({\bf Q})$ and $B({\bf Q})$ are obtained by making by substituting

[TABLE]

into $A({\bf q})$ and $B({\bf q})$ respectively, which also leads to a ‘Matsubara potential’ in Eq. (65). This transformation also diagonalizes the spring part of the ring polymer Hamiltonian in Eq. (149),

[TABLE]

Appendix C Matsubara error Liouvillian

By exploiting trigonometric identities, Eq. (67) can be given as [32]

[TABLE]

where $\hat{X}$ acts only on the non-Matsubara modes

[TABLE]

and $\hat{Y}$ acts on the Matsubara modes

[TABLE]

Although $\mathcal{L}_{\rm er}$ contains both Matsubara and non-Matsubara derivatives, expanding the trigonometric functions in Eq. (160) shows that all terms in $\mathcal{L}_{\rm er}$ contain at least one derivative in a non-Matsubara mode.

Appendix D Crossover temperature

A rough guide for the temperature beneath which quantum effects become pronounced is the crossover temperature where the first ring polymer normal mode becomes unstable, defined as [56]

[TABLE]

where $\omega_{b}$ is the imaginary frequency at the top of the barrier. Since at the maximum $dV(q)/dq=0$ by construction, the potential can be expanded as $V(q-q^{\ddagger})\simeq V(q^{\ddagger})-m\omega_{b}^{2}q^{2}/2+\mathcal{O}(q^{3})$, and $\omega_{b}$ therefore provides a guide concerning how ‘peaked’ the barrier is, as sketched in Fig. 8.