Boltzmann-conserving classical dynamics in quantum time-correlation functions: Matsubara dynamics
Timothy J. H. Hele, Michael J. Willatt, Andrea Muolo, Stuart C. Althorpe
BBoltzmann-conserving classical dynamics in quantum time-correlation functions: ‘Matsubara dynamics’
Timothy J. H. Hele, Michael J. Willatt, Andrea Muolo, a) and Stuart C. Althorpe b) Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW, UK. (Dated: 15 October 2018)
We show that a single change in the derivation of the linearized semiclassical-initial value representation(LSC-IVR or ‘classical Wigner approximation’) results in a classical dynamics which conserves the quantumBoltzmann distribution. We rederive the (standard) LSC-IVR approach by writing the (exact) quantumtime-correlation function in terms of the normal modes of a free ring-polymer (i.e. a discrete imaginary-timeFeynman path), taking the limit that the number of polymer beads N → ∞ , such that the lowest normal-modefrequencies take their ‘Matsubara’ values. The change we propose is to truncate the quantum Liouvillian, notexplicitly in powers of (cid:126) at (cid:126) (which gives back the standard LSC-IVR approximation), but in the normal-mode derivatives corresponding to the lowest Matsubara frequencies. The resulting ‘Matsubara’ dynamics is inherently classical (since all terms O ( (cid:126) ) disappear from the Matsubara Liouvillian in the limit N → ∞ ), andconserves the quantum Boltzmann distribution because the Matsubara Hamiltonian is symmetric with respectto imaginary-time translation. Numerical tests show that the Matsubara approximation to the quantum time-correlation function converges with respect to the number of modes, and gives better agreement than LSC-IVR with the exact quantum result. Matsubara dynamics is too computationally expensive to be applied tocomplex systems, but its further approximation may lead to practical methods. Copyright (2015) AmericanInstitute of Physics. This article may be downloaded for personal use only. Any other use requires priorpermission of the author and the American Institute of Physics. The following article appeared in the Journalof Chemical Physics, 142, 134103 (2015) and may be found at http://dx.doi.org/10.1063/1.4916311
I. INTRODUCTION
Dynamical properties at thermal equilibrium are ofcentral importance to chemical physics.
Sometimesthese properties can be simulated adequately by entirelyclassical means. But there are plenty of cases, e.g.the spectrum of liquid water, hydrogen-diffusion onmetals, and proton/hydride-transfer reactions, forwhich one needs to evaluate time-correlation functions ofthe form1
Z C AB ( t ) = 1 Z Tr (cid:104) e − β ˆ H ˆ Ae i ˆ Ht/ (cid:126) ˆ Be − i ˆ Ht/ (cid:126) (cid:105) (1)(where Z is the partition function, β ≡ /k B T , Tr in-dicates a complete sum over states, and the other nota-tion is defined in Sec. II). Such time-correlation functionsare already approximate, since they employ the quan-tum Boltzmann distribution e − β ˆ H /Z in place of the ex-act quantum-exchange statistics; but this approximationis usually adequate (since the thermal wavelength is typi-cally much smaller than the separations between identicalparticles). What is less well understood is the extent towhich such functions can be further approximated by re-placing the exact quantum dynamics by classical dynam-ics (whilst retaining the quantum Boltzmann statistics).The standard way to make this approximation is touse the linearized semiclassical-initial value representa-tion (LSC-IVR, sometimes called the ‘classical Wigner’ a) Current address: Lab. f¨ur Physikalische Chemie, ETH Z¨urich,CH-8093 Z¨urich, Switzerland b) Corresponding author: [email protected] approximation), in which the quantum Liouvillianis expanded as a power series in (cid:126) , then truncated at (cid:126) .Miller and later Shi and Geva showed that this ap-proximation is equivalent to linearizing the displacementbetween forward and backward Feynman paths in theexact quantum time-propagation, which removes the co-herences, thus making the dynamics classical. The LSC-IVR retains the Boltzmann quantum statistics inside aWigner transform, is exact in the zero-time, harmonicand high-temperature limits, and has been developed intoa practical method by several authors. However, ithas a serious drawback: the classical dynamics does notin general preserve the quantum Boltzmann distribution,and thus the quality of the statistics deteriorates overtime.A number of methods have been developed to getround this problem, all of which appear to some extent tobe ad hoc. Some of these methods are obtained by replac-ing the plain Newtonian dynamics in the LSC-IVR by aneffective (classical) dynamics which preserves the Boltz-mann distribution.
Others, such as the popular cen-troid molecular dynamics (CMD) and ring-polymermolecular dynamics (RPMD), are more heuris-tic (and still not fully understood) but have the advan-tage that they can be implemented directly in classi-cal molecular dynamics codes. An intriguing propertyof CMD and RPMD is that, for some model systems(e.g. the one-dimensional quartic oscillator ), thesemethods give better agreement than LSC-IVR with theexact quantum result, even though, like LSC-IVR, theycompletely neglect real-time quantum coherence.This last point suggests that the failure of LSC-IVR topreserve the quantum Boltzmann distribution may arise,not from its neglect of quantum coherence, but from its1 a r X i v : . [ phy s i c s . c h e m - ph ] A p r nclusion of ‘rogue’ components in the classical dynamics.The present paper develops a theory that supports thisspeculation. We isolate a core, Boltzmann-conserving,classical dynamics, which we call ‘Matsubara dynamics’(for reasons to be made clear). Matsubara dynamics isfar too expensive to be used as a practical method, butis likely to prove useful in understanding methods suchas CMD and RPMD, and perhaps in developing new ap-proximate methods.The paper is structured as follows. Section II gives keybackground material including the well known ‘Moyal se-ries’ derivation of the LSC-IVR. Section III re-expressesthe standard results of Sec. II in terms of ‘ring-polymer’coordinates, involving points along the imaginary-timepath-integrals that describe the quantum Boltzmannstatistics. Section IV gives the new results, showing thatsmooth Fourier-transformed combinations of the ring-polymer coordinates lead to an inherently classical dy-namics which is quantum-Boltzmann-conserving . SectionV reports numerical tests on one-dimensional models.Section VI concludes the article. II. BACKGROUND THEORY
We start by defining the terms and notation to be usedin classical and quantum Boltzmann time-correlationfunctions (IIA and IIB), and by writing out the standardMoyal-series derivation of the LSC-IVR (IIC).
A. Classical correlation functions
Without loss of generality, we can consider an F -dimensional Cartesian system with position coordinates q ≡ q , . . . , q F , momenta p , mass m and Hamiltonian H ( p , q ) = p m + V ( q ) (2)The thermal time-correlation function between observ-ables A ( p , q ), B ( p , q ) is then c AB ( t ) = 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q e − βH ( p , q ) × A ( p , q ) B ( p t , q t ) (3)where (cid:82) d p ≡ (cid:82) ∞−∞ dp . . . (cid:82) ∞−∞ dp F (and similarly for q ),and p t ≡ p t ( p , q , t ) and q t ≡ q t ( p , q , t ) are the momentaand positions after the classical dynamics has evolved fora time t .Alternatively, we can express B ( p t , q t ) as a functionof the initial phase-space coordinates ( p , q ): B ( p t , q t ) ≡ B [ p t ( p , q , t ) , q t ( p , q , t )] ≡ B ( p , q , t ) (4) such that c AB ( t ) = 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q e − βH ( p , q ) × A ( p , q ) B ( p , q , t )= 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q e − βH ( p , q ) × A ( p , q ) e L F t B ( p , q ,
0) (5)where the (classical) Liouvillian L F is L F = 1 m p · ∇ q − V ( q ) ←−∇ q · −→∇ p (6)with ∇ q = ∂∂q ... ∂∂q F (7)and the arrows indicate the direction in which the deriva-tive operator is applied (and the backward arrow indi-cates that the derivative is taken only of V ( q )—not ofany terms that may precede V ( q ) in any integral). Equa-tion (5) is less practical than Eq. (3) (which propagatesindividual trajectories rather than the distribution func-tion B ( p , q , t )) but is better for comparison with the ex-act quantum expression.An essential property of the dynamics is that it pre-serves the (classical) Boltzmann distribution, which fol-lows because H ( p , q ) is a constant of the motion. As aresult, we can rearrange Eq. (5) as c AB ( t ) = 1(2 π (cid:126) ) N (cid:90) d p (cid:90) d q e − βH ( p , q ) × (cid:2) e −L F t A ( p , q ) (cid:3) B ( p , q ,
0) (8)showing that c AB ( t ) satisfies c AB ( t ) = c BA ( − t ) (9)which is the detailed balance condition. B. Quantum correlation functions
For clarity of presentation, we will derive the results inSecs. III and IV for a one-dimensional quantum systemwith Hamiltonian ˆ H = ˆ T + ˆ V , kinetic energy operatorˆ T = ˆ p / m , potential energy operator ˆ V = V (ˆ q ), posi-tion and momentum operators ˆ q, ˆ p , and mass m . How-ever, the results we derive in Secs. III and IV are appli-cable immediately to systems with any number of dimen-sions (see Sec. IV.D).The simplest form of quantum-Boltzmann time-correlation function is that given in Eq. (1), but C AB ( t )is difficult to relate to the classical time-correlation func-tion c AB ( t ), because it does not satisfy Eq. (9) and is not2n general real. We therefore use the Kubo-transformedtime-correlation function C AB ( t ) = Tr (cid:104) K β ( ˆ A ) e i ˆ Ht/ (cid:126) ˆ Be − i ˆ Ht/ (cid:126) (cid:105) (10)with K β ( ˆ A ) = 1 β (cid:90) β dλ e − λ ˆ H ˆ Ae − ( β − λ ) ˆ H (11)This function gives an equivalent description of the dy-namics to C AB ( t ), to which it is related by a simpleFourier-transform formula. It is easy to show (by noting that e − λ ˆ H and e − i ˆ Ht/ (cid:126) commute in Eq. (10)) that C AB ( t ) satisfies the detailedbalance relation C AB ( t ) = C BA ( − t ) (12)This relation also ensures that C AB ( t ) is real (since re-versing the order of operators in the trace gives C AB ( t ) = C ∗ BA ( − t )).The t = 0 limit of C AB ( t ) can be expressed in termsof a classical Boltzmann distribution over an extendedphase space of ‘ring-polymers’. When ˆ A and ˆ B arefunctions A (ˆ q ) and B (ˆ q ) of the position operator ˆ q , thering-polymer expression is C AB (0) = lim N →∞ π (cid:126) ) N (cid:90) d p (cid:90) d q × A ( q ) B ( q ) e − β N R N ( p , q ) (13)where β N = β/N , (cid:82) d p ≡ (cid:82) ∞−∞ dp . . . (cid:82) ∞−∞ dp N and simi-larly for (cid:82) d q , and A ( q ) = 1 N N (cid:88) i =1 A ( q i ) , B ( q ) = 1 N N (cid:88) i =1 B ( q i ) (14) R N ( p , q ) = T N ( p , q ) + U N ( q ) (15) T N ( p , q ) = p m + m β N (cid:126) ) N (cid:88) i =1 ( q i +1 − q i ) (16) U N ( q ) = N (cid:88) i =1 V ( q i ) (17)Similar expressions can be obtained when ˆ A and ˆ B de-pend on the momentum operator (by inserting position-momentum Fourier-transforms). To avoid confusion, weemphasise that Eq. (13) is exact at t = 0, and that we do not assume that the ring-polymer Hamiltonian R N ( p , q )generates the dynamics at t > C. The LSC-IVR approximation
1. The Wigner-Moyal series
To derive the LSC-IVR approximation to C AB ( t ), wefollow ref. 26, expanding the exact quantum Liouvillian in powers of (cid:126) . We start by rewriting Eq. (10) as C AB ( t ) = (cid:90) ∞−∞ dq (cid:90) ∞−∞ d ∆ × (cid:104) q − ∆ / | K β ( ˆ A ) | q + ∆ / (cid:105)× (cid:104) q + ∆ / | e i ˆ Ht/ (cid:126) ˆ Be − i ˆ Ht/ (cid:126) | q − ∆ / (cid:105) (18)then insert the momentum identity δ (∆ − ∆ (cid:48) ) = 12 π (cid:126) (cid:90) ∞−∞ dp e ip (∆ − ∆ (cid:48) ) / (cid:126) (19)to obtain C AB ( t ) = 12 π (cid:126) (cid:90) ∞−∞ dq (cid:90) ∞−∞ dp × [ K β ( ˆ A )] W ( p, q ) [ ˆ B ( t )] W ( p, q ) (20)where the Wigner transforms of ˆ A and ˆ B are given by[ K β ( ˆ A )] W ( p, q ) = (cid:90) ∞−∞ d ∆ e ip ∆ / (cid:126) × (cid:104) q − ∆ / | K β ( ˆ A ) | q + ∆ / (cid:105) (21)and[ ˆ B ( t )] W ( p, q ) = (cid:90) ∞−∞ d ∆ e ip ∆ / (cid:126) × (cid:104) q − ∆ / | e i ˆ Ht/ (cid:126) ˆ Be − i ˆ Ht/ (cid:126) | q + ∆ / (cid:105) . (22)(and note that we will often suppress the ( p, q ) depen-dence of [ K β ( ˆ A )] W and [ ˆ B ( t )] W ).We then differentiate Eq. (20) with respect to t , dC AB ( t ) dt = 12 π (cid:126) (cid:90) ∞−∞ dq (cid:90) ∞−∞ dp × [ K β ( ˆ A )] W (cid:20) i (cid:126) [ ˆ H, ˆ B ( t )] (cid:21) W (23)and expand the potential-energy operator in the commu-tator in powers of ∆ to obtain (cid:20) i (cid:126) [ ˆ H, ˆ B ( t )] (cid:21) W = (cid:90) ∞−∞ d ∆ e ip ∆ / (cid:126) × ˆ (cid:96) (cid:104) q − ∆ / | ˆ B ( t ) | q + ∆ / (cid:105) (24)withˆ (cid:96) = i (cid:126) m ∂∂q ∂∂ ∆ − i (cid:126) ∞ (cid:88) λ =1 , odd λ ! ∂ λ V ( q ) ∂q λ (cid:18) ∆2 (cid:19) λ (25)Noting that each power of ∆ can be generated by anapplication of − i (cid:126) ∂/∂p , we then obtain dC AB ( t ) dt = 12 π (cid:126) (cid:90) ∞−∞ dq (cid:90) ∞−∞ dp × [ K β ( ˆ A )] W ˆ L [ ˆ B ( t )] W (26)3ithˆ L = pm ∂∂q − ∞ (cid:88) λ =1 , odd λ ! (cid:18) i (cid:126) (cid:19) λ − ∂ λ V ( q ) ∂q λ ∂ λ ∂p λ (27)This is the Moyal expansion of the quantum Liouvillian inpowers of (cid:126) . If all terms are included in the series, thenthe application of ˆ L generates the exact quantum dynam-ics (as is easily proved by working backwards through thederivation just given). A compact representation of ˆ L ,which will be useful later on is,ˆ L = pm ∂∂q − V ( q ) 2 (cid:126) sin (cid:32) ←− ∂∂q (cid:126) −→ ∂∂p (cid:33) . (28)where the arrows are defined in the same way as in Eq. (6)
2. Approximating the dynamics
To obtain the LSC-IVR one notes that Eq. (27) can bewritten ˆ L = L + O ( (cid:126) ) (29)where L is the classical Liouvillian L = pm ∂∂q − ∂V∂q ∂∂p (30)and then truncates ˆ L at (cid:126) . The LSC-IVR thus amountsto replacing the quantum dynamics by classical dynam-ics, such that C AB ( t ) is approximated by C W AB ( t ) = (cid:90) ∞−∞ dq (cid:90) ∞−∞ dp × [ K β ( ˆ A )] W ( p, q ) e L t [ ˆ B (0)] W ( p, q ) (31)or equivalently C W AB ( t ) = (cid:90) ∞−∞ dq (cid:90) ∞−∞ dp × [ K β ( ˆ A )] W ( p, q ) [ ˆ B (0)] W ( p t , q t ) (32)where ( p t , q t ) are the (classical) position and momentumat time t of a trajectory initiated at ( p, q ) at t = 0.Physical insight into the LSC-IVR is obtained by go-ing back to Eq. (25), and noting that truncating ˆ L at (cid:126) is equivalent to truncating ˆ l at ∆. Since ∆ is thedifference between the origin of a forward path that ter-minates at z (at time t ) and the terminus of a backwardpath that originates at z , it follows that truncating at∆ is equivalent to linearizing the difference between theforward and backward Feynman paths at each time-step.Hence the neglect of terms O ( (cid:126) ) is valid if the forwardand backward paths are very close together, in whichcase there are no coherence effects, and the dynamics be-comes classical. The LSC-IVR is thus exact at t = 0 (where the paths become infinitessimally short), in theharmonic limit (where the are no terms O ( (cid:126) ) in ˆ L ), andin the high temperature limit (where fluctuations in ∆efficiently dephase). Despite these positive features, LSC-IVR suffers fromthe major drawback of not preserving the quantum Boltz-mann distribution (except in one of the special limits justmentioned), since in general L [ e − β ˆ H ] W (cid:54) = 0 (33)As a result, C W AB ( t ) (cid:54) = C W BA ( − t ) (34)i.e. the LSC-IVR does not satisfy detailed balance. InSecs. III-V we will investigate why this is so. III. RING-POLYMER COORDINATES
We now recast the standard expressions of Sec. II interms of ring-polymer coordinates. No new approxima-tions are obtained, but the ring-polymer versions of theseexpressions are needed for use in Sec. IV, where they willbe used to derive the quantum-Boltzmann-conserving‘Matsubara’ dynamics.
A. Ring-polymer representation ofKubo-transformed time-correlation functions
1. Exact quantum time-correlation function
Following ref. 48 (see also refs. 17 and 41), we definethe ring-polymer quantum time-correlation function tobe C [ N ] AB ( t ) = (cid:90) d q (cid:90) d ∆ (cid:90) d z A ( q ) B ( z ) × N (cid:89) l =1 (cid:104) q l − − ∆ l − / | e − β N ˆ H | q l + ∆ l / (cid:105)× (cid:104) q l + ∆ l / | e − i ˆ Ht/ (cid:126) | z l (cid:105)× (cid:104) z l | e i ˆ Ht/ (cid:126) | q l − ∆ l / (cid:105) (35)where the functions A ( q ) and B ( z ) (with z in place of q ) are defined in Eq. (14) (and we have assumed thatˆ A and ˆ B are functions of position operators to simplifythe algebra—see Sec. IVD). It is easy to show (by notingthat N − A ( q ) and B ( z ) becomeintegrals in the limit N → ∞ ) that C AB ( t ) = lim N →∞ C [ N ] AB ( t ) (36)4 IG. 1. Schematic diagram showing the structure of the(exact) Kubo-transformed quantum time-correlation functionwhen represented in ring-polymer coordinates as in Eq. (35).The red and blue dots represent the coordinates q l and z l ;solid lines represent stretches of imaginary time of length β N (cid:126) ;arrows represent forward-backward propagations in real time. In other words, Eq. (35) in the limit N → ∞ is justan alternative way of writing out the standard Kubo-transformed time-correlation function C AB ( t ). The ad-vantage of Eq. (35) is that it emphasises the symmetry ofthe entire path-integral expression with respect to cyclicpermutations of the coordinates q l → q l +1 (see Fig. 1);this symmetry is otherwise hidden in the conventionalexpression for C AB ( t ) [Eq. (10)].
2. Ring-polymer representation of theLSC-IVR
It is straightforward to derive the LSC-IVR approx-imation from Eq. (35) by generalizing the steps inSec. IIC. We insert an identity δ (∆ l − ∆ (cid:48) l ) = 12 π (cid:126) (cid:90) ∞−∞ dp l e ip l (∆ l − ∆ (cid:48) l ) / (cid:126) (37)for each value l = 1 , . . . , N , to obtain C [ N ] AB ( t ) = 1(2 π (cid:126) ) N (cid:90) d q (cid:90) d p (cid:104) e − β ˆ H ˆ A (cid:105) N ( p , q ) (cid:104) ˆ B ( t ) (cid:105) N ( p , q ) (38)where (cid:104) e − β ˆ H ˆ A (cid:105) N ( p , q ) = (cid:90) d ∆ A ( q ) N (cid:89) l =1 (cid:104) q l − − ∆ l − / | e − β N ˆ H | q l + ∆ l / (cid:105) e ip l ∆ l / (cid:126) (39)and (cid:104) ˆ B ( t ) (cid:105) N ( p , q ) = (cid:90) d ∆ (cid:90) d z B ( z ) N (cid:89) l =1 (cid:104) q l − ∆ l / | e − i ˆ Ht/ (cid:126) | z l (cid:105)(cid:104) z l | e i ˆ Ht/ (cid:126) | q l + ∆ l / (cid:105) e ip l ∆ l / (cid:126) (40)are generalized Wigner transforms (and we will often sup-press the dependence on ( p , q ) in what follows). Notethat [ · ] N and [ · ] N have different forms: [ · ] N is a sum ofproducts of one-dimensional Wigner transforms, whereas[ · ] N is more complicated, with each product couplingvariables in l and l + 1. Note that since we have speci-fied that ˆ B is a function of just the position operator (inorder to simplify the algebra—see Sec. IVD), it followsthat (cid:104) ˆ B (0) (cid:105) N ( p , q ) = B ( q ) (41)The next step is to obtain the ring-polymer representa-tion of the (exact) quantum Liouvillian, which involvesa straightforward generalization of Eqs. (23)-(27). Wedifferentiate C [ N ] AB ( t ) with respect to t , obtain a sum of N Heisenberg time-derivatives, and expand each mem-ber in powers of ∆ l to obtain an N -fold generalization ofEqs. (24) and (25). On replacing powers of ∆ l by powers of − i (cid:126) ∂/∂p l , we obtain dC [ N ] AB ( t ) dt = 1(2 π (cid:126) ) N (cid:90) d q (cid:90) d p × (cid:104) e − β ˆ H ˆ A (cid:105) N ˆ L N (cid:104) ˆ B ( t ) (cid:105) N (42)whereˆ L N = N (cid:88) l =1 p l m ∂∂q l − V ( q l ) 2 (cid:126) sin (cid:32) ←− ∂∂q l (cid:126) −→ ∂∂p l (cid:33) . (43)and the arrow notation is as used in Eq. (6). We canwrite this expression more compactly in terms of U N ( q )in Eq. (17) asˆ L N = 1 m p · ∇ q − U N ( q ) 2 (cid:126) sin (cid:18) (cid:126) ←−∇ q · −→∇ p (cid:19) . (44)5since all mixed derivatives of U N ( q ) are zero).Following Sec. IIC, we then truncate the exact Liou-villian at (cid:126) such thatˆ L N = L N + O ( (cid:126) ) (45)with L N = N (cid:88) l =1 p l m ∂∂q l − ∂V ( q l ) ∂q l ∂∂p l (46)The ring-polymer version of LSC-IVR thus approximatesthe exact dynamics by the classical dynamics of N inde-pendent particles, each initiated at a phase-space point( p l , q l ). The ring-polymer LSC-IVR time-correlationfunction is C W[ N ] AB ( t ) = 1(2 π (cid:126) ) N (cid:90) d q (cid:90) d p × (cid:104) e − β ˆ H ˆ A (cid:105) N e L N t (cid:104) ˆ B (0) (cid:105) N = 1(2 π (cid:126) ) N (cid:90) d q (cid:90) d p × (cid:104) e − β ˆ H ˆ A (cid:105) N (cid:104) ˆ B (0) (cid:105) N ( p t , q t ) (47)where (cid:104) ˆ B (0) (cid:105) N ( p t , q t ) indicates that this Wigner trans-form takes its t = 0 form, but is expressed as a functionof the momenta and positions ( p t , q t ) of the N indepen-dent particles at time t . It is easy show (by noting thatone can integrate out N − p l ) that C W AB ( t ) = lim N →∞ C W[ N ] AB ( t ) (48)i.e. that the truncation of ˆ L N at (cid:126) gives the standardLSC-IVR approximation in the limit N → ∞ (as wouldbe expected, since we have approximated the exact quan-tum Kubo time-correlation function of Eqs. (35) and (36)by truncating the quantum Liouvillian at (cid:126) ). B. Normal mode coordinates
1. Definition
The advantage of ring-polymer coordinates is that wecan now transform to sets of global coordinates describingcollective motion of the individual coordinates ( p l , q l , ∆ l ).The choice of global coordinates is not unique. We willfind it convenient to use the normal modes of a free ring-polymer, namely the linear combinations of q l thatdiagonalize T N ( p , q ) of Eq. (16). These are simply dis-crete Fourier transforms, which for odd N (which we willassume, to simplify the algebra ), are Q n = N (cid:88) l =1 T ln q l , n = 0 , ± , . . . , ± ( N − / T ln = N − / n = 0 (cid:112) /N sin(2 πln/N ) n = 1 , . . . , ( N − / (cid:112) /N cos(2 πln/N ) n = − , . . . , − ( N − / P n in terms of p l , and D n in terms of∆ l . The associated normal frequencies take the form ω n = 2 β N (cid:126) sin (cid:16) nπN (cid:17) (51)such that the ring-polymer expression for C AB (0)[Eq. (13)] can be rewritten as C AB (0) = lim N →∞ π (cid:126) ) N (cid:90) d P (cid:90) d Q × A ( Q ) B ( Q ) e − β N R N ( P , Q ) (52)where the normal-mode expression for the ring-polymerHamiltonian R N ( P , Q ) is R N ( P , Q ) = ( N − / (cid:88) n = − ( N − / P n m + m ω n Q n + U N ( Q )(53)and A ( Q ), B ( Q ) and U N ( Q ) are obtained by making thesubstitution q l = ( N − / (cid:88) n = − ( N − / T ln Q n (54)into A ( q ), B ( q ) and U N ( q ) of Eqs. (14)-(17). Note thedefinition of the sign of ω n in Eq. (51), which results insomewhat neater expressions later on. Note also that R N ( P , Q ) will not be used to generate the dynamics inany of the expressions derived below which, like the dy-namics of Sec. IIIA, will involve N independent particlesunconnected by springs.
2. Time-correlation functions
It is straightforward to convert Eq. (38) into normalmode coordinates using the orthogonal transformationsin Eq. (54), to obtain C [ N ] AB ( t ) = 1(2 π (cid:126) ) N (cid:90) d P (cid:90) d Q × (cid:104) e − β ˆ H ˆ A (cid:105) N ( P , Q ) (cid:104) ˆ B ( t ) (cid:105) N ( P , Q ) (55)where (cid:90) d P ≡ ( N − / (cid:89) n = − ( N − / (cid:90) ∞−∞ dP n (56)6nd (cid:82) d Q is similarly defined. The generalized Wignertransforms in Eq. (55) are obtained using Eq. (54) tosubstitute ( P , Q , D ) for ( p , q , ∆ ) in Eqs. (39) and (40),and thus contain products of exp( iP n D n / (cid:126) ) in place ofexp( ip l ∆ l / (cid:126) ). At t = 0, one obtains (cid:104) ˆ B (0) (cid:105) N ( P , Q ) = B ( Q ) (57)where B ( Q ) is obtained by substituting Q for q in B ( q )of Eq. (14).The (exact) quantum dynamics is described by dC [ N ] AB ( t ) dt = 1(2 π (cid:126) ) N (cid:90) d P (cid:90) d Q × (cid:104) e − β ˆ H ˆ A (cid:105) N ˆ L N (cid:104) ˆ B ( t ) (cid:105) N (58)where the Liouvillian ˆ L N is obtained by expressing ˆ L N of Eq. (46) in terms of normal modes, which givesˆ L N = 1 m P · ∇ Q − U N ( Q ) 2 (cid:126) sin (cid:18) (cid:126) ←−∇ Q · −→∇ P (cid:19) . (59)in which U N ( Q ) is obtained by substituting Q for q in U N ( q ) of Eq. (17).As in Sec. IIIA, the LSC-IVR dynamics is obtained bytruncating ˆ L N at (cid:126) to give L N = ( N − / (cid:88) n = − ( N − / P n m ∂∂Q n − ∂U N ( Q ) ∂Q n ∂∂P n (60)after which one obtains C W[ N ] AB ( t ) in terms of normalmodes, which gives the (standard) LSC-IVR result inthe limit N → ∞ , according to Eq. (48). Hence all wehave done in Eqs. (55)-(60) is to re-express the resultsof Sec. IIIA in terms of normal mode coordinates. Theadvantages of doing this will become clear shortly. C. Matsubara modes
We now consider the M lowest frequency ring-polymernormal modes in the limit N → ∞ , such that M (cid:28) N .The frequencies ω n tend to the values (cid:101) ω n = lim N →∞ ω n = 2 nπβ (cid:126) , | n | ≤ ( M − / and so we will refer to these M modesin the limit N → ∞ as the ‘Matsubara modes’. TheMatsubara modes have the special property that any su-perposition of them produces a distribution of the coor-dinates q l which is a smooth and differentiable functionof imaginary time τ , such that q l = q ( τ ) , τ = β N (cid:126) l, l = 1 , . . . , N (62) FIG. 2. Schematic diagram showing that superpositions ofMatsubara modes give distributions of path-integral coordi-nates q l which are smooth, differentiable functions of imagi-nary time τ . Inclusion of non-Matsubara modes gives jaggeddistributions. (see Appendix A). Hence distributions made up of super-positions of the Matsubara modes resemble the sketch inFig. 2. We will often write the Matsubara modes usingthe notation (cid:101) Q n = lim N →∞ Q n √ N , n = 0 , ± , . . . , ± ( M − / (cid:101) P n , (cid:101) D n ). The extra factor of N − / ensures that (cid:101) Q n scales as N and converges in the limit N → ∞ ; e.g. (cid:101) Q is the centroid (centre of mass) ofthe smooth distribution q ( τ ). We will refer to the other N − M normal modes as the ‘non-Matsubara modes’. Ingeneral, these modes give rise to jagged (i.e. non-smooth,non-differentiable with respect to τ ) distributions of q l (see Fig. 2). Matsubara modes have a long history in path-integral descriptions of equilibrium properties, since theygive rise to an alternative ring-polymer expression for C AB (0). If we define C [ M ] AB (0) = α M π (cid:126) (cid:90) d (cid:101) P (cid:90) d (cid:101) Q A ( (cid:101) Q ) B ( (cid:101) Q ) e − β (cid:101) R M ( (cid:101) P , (cid:101) Q ) (64)7ith (cid:101) R M ( (cid:101) P , (cid:101) Q ) = ( M − / (cid:88) n = − ( M − / (cid:101) P n m + m (cid:101) ω n (cid:101) Q n + (cid:101) U M ( (cid:101) Q )(65) (cid:101) U M ( (cid:101) Q ) = lim N →∞ N N (cid:88) l =1 V ( M − / (cid:88) n = − ( M − / T ln √ N (cid:101) Q n (66) α M = (cid:126) (1 − M ) [( M − / (67)then C AB (0) = lim M →∞ M (cid:28) N C [ M ] AB (0) (68)where this limit indicates that M is allowed to tend toinfinity, subject to the condition that it is always muchsmaller than N , such that the (cid:101) Q remain Matsubaramodes. In practice, a good approximation to the exactresult is reached once (cid:101) ω ( M − / exceeds the highest fre-quency in the potential V ( q ). Equation (64) is less oftenused nowadays to compute static properties, because theconvergence with respect to M is typically slower thanthe convergence of Eq. (13) with respect to N . However, Eq. (64) tells us something interesting: TheBoltzmann factor ensures that only smooth distributionsof ( p , q ) survive in C AB ( t ) at t = 0; but at t > L N [Eq. (60)] will, in general, mix inan increasing proportion of non-smooth, non-Matsubaramodes, so that the distributions of ( p , q ) become increas-ingly jagged as time evolves. The rate at which this mix-ing occurs depends on the anharmonicity of the potential V ( q ). In the special case that V ( q ) is harmonic, there isno coupling between different normal modes, so the dis-tributions in ( p , q ) remain smooth for all time. In otherwords, smooth distributions in ( p , q ) are found in two ofthe limits (zero-time and harmonic) in which the LSC-IVR is known to be exact. IV. MATSUBARA DYNAMICSA. Definition
The results of Sec. IIIC suggest that there may be aconnection between smoothness in imaginary time andclassical dynamics. We now investigate what happens ifwe constrain an initially smooth function of phase spacecoordinates ( p , q ) to remain smooth for all (real) times t >
0. We take the (exact) quantum Liouvillian ˆ L N , andinstead of truncating at (cid:126) as in Eq. (60) (which gives theLSC-IVR), we retain all powers of (cid:126) , take the N → ∞ limit, and split ˆ L N intolim N →∞ ˆ L N = L M + lim N →∞ ˆ L error ( N, M ) (69) where the ‘Matsubara Liouvillian’ L M = lim N →∞ ( M − / (cid:88) n = − ( M − / P n m ∂∂Q n − U N ( Q ) 2 (cid:126) sin ( M − / (cid:88) n = − ( M − / (cid:126) ←− ∂∂Q n −→ ∂∂P n (70)contains all terms in which the derivatives involve only the Matsubara modes, and ˆ L error ( N, M ) contains the restof the terms (given in Appendix B). We then discardˆ L error ( N, M ), approximating ˆ L N by L M . We will re-fer to the (approximate) dynamics generated by L M as‘Matsubara dynamics’. By construction, Matsubara dy-namics ensures that a distribution of ( p , q ) which is asmooth and differentiable function of τ at t = 0 will re-main so for all t > C [ M ] AB ( t ) = lim N →∞ π (cid:126) ) N (cid:90) d P (cid:90) d Q × (cid:104) e − β ˆ H ˆ A (cid:105) N e L M t (cid:104) ˆ B (0) (cid:105) N (71)We can obtain an explicit form for C [ M ] AB ( t ) by taking thesame limit as in Eq. (68), allowing M to tend to infinity,subject to M (cid:28) N , which gives (see Appendix C) C Mats AB ( t ) = lim M →∞ M (cid:28) N C [ M ] AB ( t ) (72)where C [ M ] AB ( t ) = α M π (cid:126) (cid:90) d (cid:101) P (cid:90) d (cid:101) Q A ( (cid:101) Q ) e − β [ (cid:101) H M ( (cid:101) P , (cid:101) Q ) − iθ M ( (cid:101) P , (cid:101) Q )] × e L M t B ( (cid:101) Q ) (73)in which the Matsubara Hamiltonian is (cid:101) H M ( (cid:101) P , (cid:101) Q ) = (cid:101) P m + (cid:101) U M ( (cid:101) Q ) (74)and the phase factor is θ M ( (cid:101) P , (cid:101) Q ) = ( M − / (cid:88) n = − ( M − / (cid:101) P n (cid:101) ω n (cid:101) Q − n (75)with α M , (cid:101) ω n , (cid:101) P and (cid:101) Q defined in Sec. IIIC. Note that, inderiving these equations (in Appendix C), we have notproved that C [ M ] AB ( t ) converges with M for t > M ). Wetest this convergence numerically in Sec. V.Thus when the exact dynamics is approximated byMatsubara dynamics, the quantum Boltzmann distribu-tion takes the simple form of a classical Boltzmann dis-tribution multiplied by a phase factor. At t = 0, one8ay analytically continue the phase factor (by making P n → P n − imω n Q − n ) to recover the ring-polymer dis-tribution in Eq. (64). However, it is not known whetherthis analytic continuation is valid at t > B. Matsubara dynamics is classical
We now rewrite L M in terms of ( (cid:101) P , (cid:101) Q ), to make ex-plicit its dependence on N , and we also assume that M is sufficiently large that Eq. (73) holds, allowing us toreplace U N ( Q ) /N by (cid:101) U M ( (cid:101) Q ). This gives L M = lim N →∞ m (cid:101) P · ∇ (cid:101) Q − (cid:101) U M ( (cid:101) Q ) 2 N (cid:126) sin (cid:18) (cid:126) N ←−∇ (cid:101) Q · −→∇ (cid:101) P (cid:19) . (76)In other words, the Moyal series in Matsubara space isan expansion in terms of ( (cid:126) /N ) , rather than (cid:126) . Now,it is well known that the smallness of (cid:126) cannot in gen-eral be used to justify truncating the (standard LSC-IVR) Moyal series of Eq. (27) at (cid:126) , since at least one ofthe Wigner transforms in the time-correlation function[Eq. (20)] contains derivatives that scale as (cid:126) − . How-ever, it is easy to show that the derivatives of all termsin the integral in Eq. (73) scale as N . As a result, itfollows that all derivatives higher than first order in L M vanish in the limit N → ∞ , with the result that L M = ( M − / (cid:88) n = − ( M − / (cid:101) P n m ∂∂ (cid:101) Q n − ∂ (cid:101) U M ( (cid:101) Q ) ∂ (cid:101) Q n ∂∂ (cid:101) P n (77)In other words, Matsubara dynamics is classical .This is a surprising result, which needs to be inter-preted with caution. It does not mean that the depen-dence of B ( Q ) on the Matsubara modes evolves classi-cally in the exact quantum dynamics, since the exactLiouvillian ˆ L N contains derivative terms that couple theMatsubara modes with the non-Matsubara modes (forwhich the higher-order derivatives cannot be neglected):it means that the dynamics of the Matsubara modes be-comes classical when they are decoupled from the non-Matsubara modes.One way to understand the origin of the (cid:126) /N inEq. (76) is to note that the Fourier transform between (cid:101) P n and (cid:101) D n (in the Wigner transforms of Eqs. (39) and(40)) is exp( iN (cid:101) P n (cid:101) D n / (cid:126) ). Hence the effective Planck’sconstant associated with motion in the Matsubara coor-dinates tends to zero in the limit N → ∞ . Note thatthe dependence of the Boltzmann distribution on thenon-Matsubara modes is more complicated than that ofEq. (73), and contains powers of ( (cid:126) /N ) − which cancelout the powers of ( (cid:126) /N ) in ˆ L N (which must obviously happen, since we know that the exact dynamics is not ingeneral classical).Matsubara dynamics thus has many features in com-mon with LSC-IVR: it is exact in the t = 0 limit (whenall distributions of ( p , q ) are smooth superpositions ofMatsubara modes), in the harmonic limit (where the dy-namics of the Matsubara modes is decoupled from thatof the non-Matsubara modes), and in the classical limit(since setting M = 0 in Eq. (73) gives the classical time-correlation function); and it neglects all terms O ( (cid:126) )in the (exact) quantum Liouvillian. However, Matsub-ara dynamics differs from LSC-IVR in that it also ne-glects the terms O ( (cid:126) ) that contain derivatives in thenon-Matsubara modes. One can thus regard Matsubaradynamics as a filtered version of LSC-IVR, in which theparts of the dynamics that cause the smooth distribu-tions of ( p , q ) to become jagged have been removed. C. Conservation of the quantum Boltzmanndistribution
Confining the dynamics to the space of Matsubaramodes has a major effect on the symmetry of the Hamil-tonian. The LSC-IVR Hamiltonian H N ( P , Q ) is simplythe classical Hamiltonian of N independent particles, andis thus symmetric with respect to any permutation of thephase space coordinates [e.g. ( p , q ) ↔ ( p , q )]. On re-stricting the dynamics to the Matsubara modes, mostof these symmetries are lost (since individual permuta-tions would destroy the smoothness of the distributionsof ( p , q )). However, one operation which is retained is symmetry with respect to cyclic permutation of thecoordinates, which, on restricting the dynamics to Mat-subara space, becomes a continuous, differentiable sym-metry, namely invariance with respect to translation inimaginary time: d (cid:101) H M ( (cid:101) P , (cid:101) Q ) dτ = 0 (78)(see Appendix A). It thus follows from Noether’stheorem, that d (cid:101) Λ M ( (cid:101) P , (cid:101) Q ) dτ = ddt ( M − / (cid:88) n = − ( M − / (cid:101) P n d (cid:101) Q n dτ = 0 (79)where Λ M ( (cid:101) P , (cid:101) Q ) is the Matsubara Lagrangian. In otherwords, in Matsubara dynamics, there exists a constantof the motion (in addition to the total energy) which isgiven by the term in brackets above.In Appendix A, it is shown that the phase θ M ( (cid:101) P , (cid:101) Q )in the quantum Boltzmann distribution [Eqs. (73)-(75)]can be written θ M ( (cid:101) P , (cid:101) Q ) = − ( M − / (cid:88) n = − ( M − / (cid:101) P n d (cid:101) Q n dτ (80)9nd is thus the constant of the motion associated with theinvariance of (cid:101) H M ( (cid:101) P , (cid:101) Q ) to imaginary time-translation.Since (cid:101) H M ( (cid:101) P , (cid:101) Q ) is of course also a constant of the mo-tion, it follows that Matsubara dynamics conserves thequantum Boltzmann distribution .As a result, Matsubara dynamics satisfies the detailedbalance relation C [ M ] AB ( t ) = C [ M ] BA ( − t ) (81)and gives expectation values (cid:10) ˆ B (cid:11) [ M ] ( t ) = α M π (cid:126) (cid:90) d (cid:101) P (cid:90) d (cid:101) Q × e − β [ (cid:101) H M ( (cid:101) P , (cid:101) Q ) − iθ M ( (cid:101) P , (cid:101) Q )] B ( (cid:101) Q t )= α M π (cid:126) (cid:90) d (cid:101) P t (cid:90) d (cid:101) Q t × e − β [ (cid:101) H M ( (cid:101) P t , (cid:101) Q t ) − iθ M ( (cid:101) P t , (cid:101) Q t )] B ( (cid:101) Q t )= α M π (cid:126) (cid:90) d (cid:101) P t (cid:90) d (cid:101) Q t × e − β (cid:101) R M ( (cid:101) P t , (cid:101) Q t ) B ( (cid:101) Q t )= (cid:10) ˆ B (cid:11) [ M ] (0) (82)which are independent of time (and equal to the exactquantum result in the limit M → ∞ ; see Eq. (68)). Notethat the step between the second and third lines followsfrom analytic continuation ( P n → P n − imω n Q − n ).We thus have the surprising result that a purely clas-sical dynamics (Matsubara dynamics) which uses thesmoothed Hamiltonian that arises naturally when thespace is restricted to Matsubara modes, conserves thequantum Boltzmann distribution. At first sight this mayappear counter-intuitive. For example, it is clear that theclassical dynamics will not respect zero-point energy con-straints, nor will it be capable of tunnelling. However, itis the phase θ M ( (cid:101) P , (cid:101) Q ) which converts what would other-wise be a classical Boltzmann distribution in an extendedphase-space into a quantum Boltzmann distribution, andthe phase is conserved. D. Generalizations
The derivations above can easily be generalized to sys-tems with any number of dimensions. For a system whoseclassical Hamiltonian resembles Eq. (2), there are F × M Matsubara modes, one set of M modes in each dimension.All the steps in Secs. III and IV.A-C are then the same,except that, with F dimensions instead of one, there isnow a sum of F phase terms, each resembling θ M ( (cid:101) P , (cid:101) Q ).Noether’s theorem shows that the sum of these terms andhence the quantum Boltzmann distribution is conserved.We emphasise that the derivations above were carriedout for operators ˆ A and ˆ B in C AB ( t ) which are general functions of the coordinate operators ˆ q . Matsubara dy-namics is therefore not limited to correlation functions -0.4-0.200.20.40.6 0 5 10 15 20 C [ M ] qq ( t ) t (a.u.) FIG. 3. Convergence with respect to number of modes M of the Matsubara position auto-correlation function C [ M ] qq ( t ),calculated for the quartic potential of Eq. (83), at a reciprocaltemperature of β = 2 a.u. The red lines correspond to M = 1(dots), 3 (chains), 5 (dashes) and 7 (solid). The solid blackline is the exact quantum result. -8-4048 0 5 10 15 20 θ M (cid:16) f P , f Q (cid:17) t (a.u.) FIG. 4. Evolution of the phase θ M ( (cid:101) P , (cid:101) Q ) along a single clas-sical trajectory on the quartic potential, with M = 5, and N = 5 (dots), 9 (dashes) and ∞ (solid line). The lattercorresponds to Matsubara dynamics in which the phase isconserved. involving linear operators of position. The derivationscan also be repeated, with minor modifications in the al-gebra, for the case that ˆ A and ˆ B are general functions ofthe momentum operator (which results in functions of (cid:101) P appearing in the generalised Wigner transforms). V. NUMERICAL TESTS OF THEEFFICACY OF MATSUBARA DYNAMICS
So far we have made no attempt to justify the use ofMatsubara dynamics, beyond pointing out that it is exactin all the limits in which LSC-IVR is exact, but that, un-like LSC-IVR, it also conserves the quantum Boltzmanndistribution. Here we investigate whether Matsubara dy-namics converges with respect to the number of modes M , and make numerical comparisons with the LSC-IVR,CMD and RPMD methods.10 .550.650.750.85 0 5 10 15 20 25 30 D q E t (a.u.)0.450.550.65 0 5 10 15 20 D q E t (a.u.)b)a) FIG. 5. Time-dependence of the thermal expectation value (cid:10) q (cid:11) ( t ) for the quartic potential at β = 2, computed us-ing LSC-IVR (blue), and Matsubara dynamics (red: M = 1(dots), 3 (dashes), 5 (solid)), and compared with the exactquantum result (black). The presence of the phase θ M ( (cid:101) P , (cid:101) Q ) in the Boltzmanndistribution [Eq. (73)] means that Matsubara dynamicssuffers from the sign problem, and thus cannot be usedas a practical method. However, we were able to evaluate C [ M ] qq ( t ) (i.e. C [ M ] AB ( t ) of Eq. (73) with ˆ A = ˆ q , ˆ B = ˆ q ) forsome one-dimensional model systems. For consistencywith previous work, we considered the quartic po-tential V ( q ) = 14 q (83)and the weakly anharmonic potential V ( q ) = 12 q + 110 q + 1100 q (84)where atomic units are used with m = 1. Calculationsusing potentials with intermediate levels of anharmonic-ity were found to give similar results (and are not shownhere).Figure 3 shows C [ M ] qq ( t ) for the quartic potential, at aninverse temperature of β = 2 a.u., for various values of M . These results were obtained by propagating classi-cal trajectories using the Matsubara potential (cid:101) U M ( (cid:101) Q ) togenerate the forces, subject to the Anderson thermostat (according to which each (cid:101) P n was reassigned to a valuedrawn at random from the classical Boltzmann distribu-tion every 2 atomic time units); (cid:101) U M ( (cid:101) Q ) was computed -0.6-0.300.30.6 0 5 10 15 20 25 30 C [ M ] qq ( t ) t (a.u.)-0.4-0.200.20.4 0 5 10 15 20 C [ M ] qq ( t ) t (a.u.) QuantumLSC-IVRCMDRPMDMatsubarab) QuantumLSC-IVRCMDRPMDMatsubaraa) QuantumLSC-IVRCMDRPMDMatsubara FIG. 6. Comparisons of position-autocorrelation functionscomputed using different levels of theory, for (a) the quarticpotential and (b) the weakly anharmonic potential of Eq. (84).The Matsubara results were obtained using M = 7 (quartic)and M = 5 (weakly anharmonic). by taking the N → ∞ limit analytically, as described inthe supplemental material. A total of 10 Monte Carlopoints was found necessary to converge C [ M ] qq ( t ). Extend-ing these calculations beyond M = 7 was prohibitivelyexpensive, and the final few M were particularly diffi-cult to converge (since θ M ( (cid:101) P , (cid:101) Q ) becomes increasinglyoscillatory as (cid:101) ω n increases). Nevertheless, the results inFig. 3 are sufficient to show that C [ M ] qq ( t ) converges withrespect to M , although the convergence appears to be-come slower as t increases. For the weakly anharmonicpotential, convergence to within graphical accuracy wasobtained using M = 5 for β = 2 a.u.We also confirmed numerically that Matsubara dynam-ics conserves the quantum Boltzmann distribution. Fig-ure 4 shows the phase θ M ( (cid:101) P , (cid:101) Q ) as a function of timealong a Matsubara trajectory. When a coarse number ofpolymer beads ( N = 5) is used, such that the M lowest-frequency modes are a poor approximation to the Mat-subara modes, the phase is not conserved; however, as N is increased, the variation of the phase along the tra-jectory flattens, becoming completely time-independentin the limit N → ∞ . Figure 5 plots the expectationvalue (cid:10) q (cid:11) [ M ] ( t ), which is found to be time-independentas expected from Eq. (82).Figure 6 compares the Matsubara correlation functions C [ M ] qq ( t ) for both potentials with exact quantum, LSC-11VR, CMD and RPMD results. The quartic potential at β = 2 (panel a) is a severe test for which any methodthat neglects real-time coherence fails after a single re-currence. Nevertheless, we see that Matsubara dynamicsgives a much better treatment than LSC-IVR, reproduc-ing almost perfectly the first recurrence at 6 a.u., anddamping to zero more slowly. The Matsubara resultis also better than both the CMD and RPMD results.The same trends are found for the weakly anharmonicpotential (Fig. 6, panel b), and were also found for thepotentials with intermediate anharmonicity (not shown).
VI. CONCLUSIONS
We have found that a single change in the derivationof LSC-IVR dynamics gives rise to a classical dynam-ics (‘Matsubara dynamics’) which preserves the quantumBoltzmann distribution. This change involves no explicittruncation in powers of (cid:126) , but instead a decoupling ofa subspace of ring-polymer normal modes (the Matsub-ara modes) from the other modes. The dynamics in thisrestricted space is found to be purely classical and to en-sure that smooth distributions of phase-space points (asa function of imaginary time), which are present in theBoltzmann distribution at time t = 0, remain smooth atall later times. The LSC-IVR dynamics, by contrast, in-cludes all the modes, which has the effect of breaking upthese smooth distributions, and thus failing to preservethe quantum Boltzmann distribution. Numerical testsshow that Matsubara dynamics gives consistently betteragreement than LSC-IVR with the exact quantum time-correlation functions.These results suggest that Matsubara dynamics isa better way than LSC-IVR, at least in principle, toaccount for the classical mechanics in quantum time-correlation functions. We suspect that Matsubaradynamics may be equivalent to expanding the time-dependence of the quantum time-correlation function inpowers of (cid:126) and truncating it at (cid:126) ; this is in con-trast to LSC-IVR, in which one truncates the quantumLiouvillian at (cid:126) . However, further work will be neededto prove or disprove this conjecture.Matsubara dynamics is far too expensive to be usefulas a practical method. However, it is probably a goodstarting point from which to make further approxima-tions in order to develop such methods. The numericaltests reported here show that Matsubara dynamics givesconsistently better results than both CMD and RPMD,suggesting that these popular methods may be approxi-mations to Matsubara dynamics. ACKNOWLEDGMENTS
TJHH, MJW and SCA acknowledge funding from theUK Science and Engineering Research Council. AM ac-knowledges the European Lifelong Learning Programme (LLP) for an Erasmus student placement scholarship.TJHH also acknowledges a Research Fellowship from Je-sus College, Cambridge and helpful discussions with DrAdam Harper.
Appendix A: Differentiability with respect toimaginary time
A distribution of ring-polymer coordinates q l , l =1 , . . . , N , can be written as a smooth and differentiablefunction of the imaginary time τ (0 ≤ τ < β (cid:126) ) if thelimit dq ( τ ) dτ = lim N →∞ q l +1 − q l − β N (cid:126) , τ = l (cid:126) β N (A1)exists, i.e. if lim N →∞ q l +1 − q l − ∼ N − (A2)For a distribution formed by superposing only the Mat-subara modes, we can use trigonometric identities andthe definitions in Sec. III to write q l +1 − q l − = 2 √ ( M − / (cid:88) n =1 (cid:20) cos (cid:18) πnlN (cid:19) (cid:101) Q n − sin (cid:18) πnlN (cid:19) (cid:101) Q − n (cid:21) × sin (cid:18) πnN (cid:19) (A3)Since n (cid:28) N , the sine function on the right ensures thatEq. (A2) is satisfied; also, repetition of this procedureshows that higher-order differences of order λ scale as N − λ . Hence a distribution q l formed from a superposi-tion of Matsubara modes is a smooth and differentiablefunction of τ . The same is true for distributions in p l and∆ l .To prove that the Matsubara Hamiltonian is invariantunder imaginary-time translation [Eq. (78)], we first dif-ferentiate the Matsubara potential U M ( (cid:101) Q ) with respectto τ , which gives d (cid:101) U M ( (cid:101) Q ) dτ = lim N →∞ (cid:101) P (cid:101) U M ( q ) − (cid:101) U M ( q ) β N (cid:126) (A4)where (cid:101) U M ( q ) = N (cid:88) l =1 V N (cid:88) m =1 ( M − / (cid:88) n = − ( M − / T ln T mn q m (A5)and P represents a cyclic permutation of the coordinates q m → q m +1 , such that P (cid:101) U M ( q ) = N (cid:88) l =1 V N (cid:88) m =1 ( M − / (cid:88) n = − ( M − / T ln T ( m − n q m (A6)12e then rearrange the sum over n in Eq. (A6) into T l T ( m −
1) 0 + ( M − / (cid:88) n =1 (cid:2) T ln T ( m − n + T l − n T ( m − − n (cid:3) (A7)and use trigonometric identities to show that T ln T ( m − n + T l − n T ( m − − n = T ( l +1) n T mn + T ( l +1) − n T m − n (A8)Re-ordering the sum over l and using the property that T l = N − / gives P (cid:101) U M ( q ) = (cid:101) U M ( q ) (A9)which proves that d (cid:101) U M ( (cid:101) Q ) dτ = 0 (A10)The same line of argument can be applied to the kineticenergy (cid:101) P / m , thus proving Eq. (78).To obtain the derivative of (cid:101) Q n with respect to τ (needed to prove Eq. (80)), we write d (cid:101) Q n dτ = lim N →∞ √ N N (cid:88) l =1 T ln q l +1 − q l − β N (cid:126) = lim N →∞ √ N N (cid:88) l =1 [ T ( l − n − T ( l +1) n ] q l β N (cid:126) (A11)and use trigonometric identities to obtain T l +1 n − T l − n = 2 T l − n sin(2 nπ/N ) (A12)Since n (cid:28) N , it follows that d (cid:101) Q n dτ = − (cid:101) ω n (cid:101) Q − n (A13)where (cid:101) ω n is the Matsubara frequency defined in Eq. (61). Appendix B: Error term for MatsubaraLiouvillian
The error term ˆ L error ( N, M ) of Eq. (69) is the differ-ence ˆ L N − L M between the exact quantum Liouvillianand the Matsubara Liouvillian. Using Eqs. (59) and (70)and the trigonometric identitysin( a + b ) − sin a ≡ (cid:18) b (cid:19) cos (cid:18) a + b (cid:19) (B1)we can writeˆ L error ( N, M ) = ( N − / (cid:88) n =( M +1) / P − n m ∂∂Q − n + P n m ∂∂Q n − (cid:126) U ( Q ) sin (cid:32) ˆ X (cid:33) cos (cid:32) ˆ Y + ˆ X (cid:33) (B2) with ˆ X = (cid:126) ( N − / (cid:88) n =( M +1) / ←− ∂∂Q − n −→ ∂∂P − n + ←− ∂∂Q n −→ ∂∂P n (B3)and ˆ Y = (cid:126) ( M − / (cid:88) n = − ( M − / ←− ∂∂Q n −→ ∂∂P n (B4) Appendix C: Derivation of Matsubaratime-correlation function
To obtain the expression for C [ M ] AB ( t ) in Eq. (73),we note that B ( P , Q , t ) is independent of the non-Matsubara P modes (since, by construction, these modesare not involved in the Matsubara dynamics) which cantherefore be integrated out, giving a product of Diracdelta-functions in the non-Matsubara (cid:101) D modes. As aresult, the Wigner transform (cid:104) e − β ˆ H ˆ A (cid:105) N in Eq. (71) re-duces to (cid:104) e − β ˆ H ˆ A (cid:105) N ( P M , Q )= (2 π (cid:126) ) N − M A ( Q ) (cid:90) d D M ( M − / (cid:89) n = − ( M − / e iP n D n / (cid:126) × N (cid:89) l =1 (cid:104) η − l − ( Q , D M ) | e − β N ˆ H | η + l ( Q , D M ) (cid:105) (C1)where P M and D M include only the Matsubara modes(and Q includes all N modes), and η ± l ( Q , D M ) = ( N − / (cid:88) n = − ( N − / T ln Q n ± ( M − / (cid:88) n = − ( M − / T ln D n / η ± l on ( Q , D M ) will be sup-pressed in what follows). Expressing the bra-ket in ring-polymer form, and using trigonometric identities, we ob-tain (cid:104) e − β ˆ H ˆ A (cid:105) N ( P M , Q )= (2 π (cid:126) ) N − M (cid:18) m πβ N (cid:126) (cid:19) N/ A ( Q ) (cid:90) d D M × e − β N mf M ( Q , D M ) / M − / (cid:89) n = − ( M − / e iP n D n / (cid:126) × exp (cid:40) − β N (cid:34) N (cid:88) l =1 V ( η − l ) + V ( η + l ) (cid:35)(cid:41) (C3)13here f M ( Q , D M )= 4( β N (cid:126) ) M − / (cid:88) n = − ( M − / (cid:18) Q n sin nπN + D − n nπN (cid:19) + ( N − / (cid:88) n =( M +1) / ( Q n + Q − n ) ω n (C4)On taking the limit N → ∞ , and converting D M to (cid:101) D , we find that the Gaussians involving D M in Eq. (C3)have the form exp (cid:16) − m (cid:101) D n N / β (cid:126) (cid:17) (C5)i.e. each Gaussian in (cid:101) D becomes a Dirac delta-functionin the limit N → ∞ . This allows us to replace the thirdline in Eq. (C3) byexp − β N N (cid:88) l =1 V ( N − / (cid:88) n = − ( N − / T ln Q n (C6)and to integrate out the (cid:101) D , giving (cid:104) e − β ˆ H ˆ A (cid:105) N ( P M , Q )= (cid:18) πmβ N (cid:19) ( N − M ) / A ( Q ) × e − β N P M / m ( M − / (cid:89) n = − ( M − / e iP n Q − n tan( nπ/N ) / (cid:126) × exp − β N N (cid:88) l =1 V ( N − / (cid:88) n = − ( N − / T ln Q n × exp − β N m ( N − / (cid:88) n =( M +1) / ( Q n + Q − n ) ω n (C7)We then substitute this expression into the integral ofEq. (71) (with (cid:82) d P replaced by (cid:82) d P M ), and take thelimit M → ∞ (subject to M (cid:28) N ), which allows us tointegrate out the non-Matsubara modes in Q . Use of theformula N − (cid:89) n =1 sin ( nπ/N ) = N/ N − (C8)then gives Eqs. (73)-(75). REFERENCES D. Chandler,
Introduction to Modern Statistical Me-chanics (Oxford University Press, New York, 1987). D. Frenkel and B. Smit,
Understanding Molecular Sim-ulation (Academic Press, London, 2002). J. Liu, W.H. Miller, F. Paesani, W. Zhang and D.A.Case, J. Chem. Phys. , 164509 (2009). S.D. Ivanov, A Witt, M. Shiga and D. Marx, J. Chem.Phys. , 031101 (2010). S. Habershon, G.S. Fanourgakis and D.E. Manolopou-los, J. Chem. Phys. , 074501 (2008). E.M. McIntosh, K.T. Wikfeldt, J. Ellis, A. Michaelidesand W. Allison, J. Phys. Chem. Lett. Y.V. Suleimanov, J. Phys. Chem. A , 11141 (2012). N. Boekelheide, R. Salom´on-Ferrer and T.F. Miller III,PNAS , 16159 (2011). R. Collepardo-Guevara, I.R. Craig andD.E. Manolopoulos, J. Chem. Phys. , 144502(2008). S. Hammes-Schiffer and A.A. Stuchebrukhov, Chem.Rev. , 6939 (2010). M. Topaler and N. Makri, J. Chem. Phys. , 7500(1994). Q. Shi, L. Zhu and L. Chen, J. Chem. Phys. ,044505 (2011). J.B. Rommel, T.P.M. Goumans and J. K¨astner, J.Chem. Theor. Comput. , 690 (2011). This definition of the time-correlation function (exclud-ing the 1 /Z ) will be used throughout. W.H. Miller, J. Phys. Chem. A , 2942 (2001). H. Wang, X. Sun and W.H. Miller, J. Chem. Phys. ,9726 (1998). Q. Shi and E. Geva, J. Chem. Phys. , 8173 (2003). T. Yamamoto, H. Wang and W.H. Miller, J. Chem.Phys. , 7335 (2002). J. Liu and W.H. Miller, J. Chem. Phys. , 074113(2009). J.A. Poulsen, G. Nyman and P.J. Rossky, J. Chem.Phys. , 12179 (2003). Q. Shi and E. Geva, J. Phys. Chem. A , 9059(2003). J. Liu, “Recent advances in the linearized semiclassi-cal initial value representation/classical Wigner modelfor the thermal correlation function,” Int. J. Quan-tum Chem. (published online). J. Beutier, D. Borgis, R. Vuilleumier and S. Bonella, J.Chem. Phys. , 084102 (2014). S. Bonella and D.F. Coker, J. Chem. Phys. , 194102(2005). P. Huo, T.F. Miller III and D.F. Coker, J. Chem. Phys. , 151103 (2013). M. Hillery, R.F. O’Connell, M.O. Scully andE.P. Wigner, Phys. Rep. , 121 (1984). E.J. Heller, J. Chem. Phys. , 1289 (1976). J. Liu and W.H. Miller, J. Chem. Phys. , 104102(2011). J. Liu, J. Chem. Phys. , 224107 (2014). J.A. Poulsen, personal communication (2014). S. Jang and G.A. Voth, J. Chem. Phys. , 2371(1999). T.D. Hone, P.J. Rossky and G.A. Voth, J. Chem. Phys.14 , 154103 (2006). I.R. Craig and D.E. Manolopoulos, J. Chem. Phys. , 3368 (2004). I.R. Craig and D.E. Manolopoulos, J. Chem. Phys. , 084106 (2005). I.R. Craig and D.E. Manolopoulos, J. Chem. Phys. , 034102 (2005). T.E. Markland and D.E. Manolopoulos, J. Chem. Phys. , 024105 (2008). S. Habershon, D.E. Manolopoulos, T.E. Markland andT.F. Miller III, Annu. Rev. Phys. Chem. , 387(2013). A.R. Menzeleev, N. Ananth and T.F. Miller III, J.Chem. Phys. , 074106 (2011). J.S. Kretchmer and T.F. Miller III, J. Chem. Phys. , 134109 (2013). A.R. Menzeleev, F. Bell and T.F. Miller III, J. Chem.Phys. , 064103 (2014). N. Ananth, J. Chem. Phys. , 124102 (2013). Y. Li, Y.V. Suleimanov, M. Yang, W.H. Green andH. Guo, J. Phys. Chem. Lett. , 48 (2013). R. P´erez de Tudela, F.J. Aoiz, Y.V. Suleimanov andD.E. Manolopoulos, J. Phys. Chem. Lett. , 493 (2012). Y.V. Suleimanov, W.J. Kong, H. Guo and W.H. Green, , 244103 (2014). P.E. Videla, P.J. Rossky and D. Laria, J. Chem. Phys. , 174315 (2013). M. Rossi, M. Ceriotti and D.E. Manolopoulos, J.Chem. Phys. , 234116 (2014). J.O. Richardson and S.C. Althorpe, J. Chem. Phys. , 214106 (2009). T.J.H. Hele and S.C. Althorpe, J. Chem. Phys. ,084108 (2013). S.C. Althorpe and T.J.H. Hele, J. Chem. Phys. ,084115 (2013). T.J.H. Hele and S.C. Althorpe, J. Chem. Phys. ,084116 (2013). T.J.H. Hele,
Quantum Transition-State Theory , PhDThesis (University of Cambridge, 2014). J.O. Richardson and M. Thoss, J. Chem. Phys. ,031102 (2013). Y. Zhang, T. Stecher, M.T. Cvitas and S.C. Althorpe,J. Phys. Chem. Lett. , 3976 (2014). R. Zwanzig,
Nonequilibrium Statistical Mechanics ,(Oxford University Press, New York, 2001). D. Chandler and P.G. Wolynes, J. Chem. Phys. ,4078 (1981). M. Parrinello and A. Rahman, J. Chem. Phys. , 860(1984). D.M. Ceperley, Rev. Mod. Phys. , 279 (1995). C. Chakravarty, Int. Rev. Phys. Chem. , 421 (1997). This paragraph is a heuristic summary of what is de-rived properly in refs. 16 and 17. One could alternatively define the coordinates q and ∆ such that [ · ] N has the form of [ · ] N and vice versa; this would yield identical results in the limit N → ∞ . All the derivations reported here can also be done foreven N , at the cost of doubling the amount of algebra inorder to deal with the awkward N/ T. Matsubara, Prog. Theor. Phys. , 351 (1955). Clearly this distinction is artificial, since, for any M (cid:28) N , there will be ‘non-Matsubara’ modes for which | n | > ( M − / | n | (cid:28) N , and which therefore also satisfyEq. (61). However, all we need to know is that all ofthe M Matsubara modes become increasingly smoothin the limit N → ∞ , and that the majority of the N − M non-Matsubara modes do not. D.L. Freeman and J.D. Doll, J. Chem. Phys. , 5709(1984). C. Chakravarty, M.C. Gordillo and D.M. Ceperley, J.Chem. Phys. , 2123 (1998). The (cid:126) /N -dependence in Eq. (76) is not an artifact ofhaving scaled ( P , Q ) to ( (cid:101) P , (cid:101) Q ); each derivative with re-spect to P n or Q n in Eq. (70) carries an implicit scalingof N − / . We could thus have derived Matsubara dynamics bystarting with the LSC-IVR Liouvillian L N of Eq. (60),then discarding the non-Matsubara derivative terms;but this would have hidden the important propertythat Matsubara dynamics is inherently classical. Reflection symmetries (e.g. τ → β (cid:126) − τ ) are also re-tained. H. Goldstein,
Classical Mechanics , 2nd ed. (Addison-Wesley, Reading, Massachusetts, 1980). See supplemental material at [URL] for further detailsof these numerical calculations. The persistence of the oscillations in the Matsubararesult for the quartic oscillator (Fig. 6a) suggests thata Matsubara version of (non-linearized) semiclassical-IVR (i.e. the next approximation up in the hierarchy,in which individual forward-backward paths are treatedclassically and assigned phases ) may give very closeagreement with the exact quantum result for signifi-cantly longer times. These two approximations are identical only in the lim-its in which LSC-IVR and Matsubara dynamics agreeand are both exact (i.e. the short-time, harmonic andhigh-temperature limits). Note that if ˆ B is a function of the momentum opera-tor, then B ( P , Q , t ) does depend on the non-Matsubara P coordinates, but that this dependence has a known,analytic, form, such that these coordinates can be inte-grated out, converting the original dependence on thenon-Matsubara P coordinates into a dependence on thenon-Matsubara Q coordinates. I.S. Gradshteyn and I.M. Ryzhik,