Path-integral approximations to quantum dynamics

Imaginary-time path-integral or ‘ring-polymer’ methods have been used to simulate quantum (Boltzmann) statistical properties since the 1980s. This article reviews the more recent extension of such methods to simulate quantum dynamics, summarising the chain of approximations that links practical path-integral methods, such as centroid molecular dynamics (CMD) and ring-polymer molecular dynamics (RPMD), to the exact quantum Kubo time-correlation function. We focus on single-surface Born–Oppenheimer dynamics, using the infrared spectrum of water as an illustrative example, but also survey other recent applications and practical techniques, as well as the limitations of current methods and their scope for future development.


Introduction
Since the 1980s, imaginary-time path-integrals [1] or 'ring-polymers' have been exploited as a practical technique for simulating quantum Boltzmann statistical properties [2][3][4][5] (and in some cases Bosonic statistical properties [5][6][7]). 1 This review addresses the more recent extension of such methods to calculate the dynamical properties of atomic nuclei in systems which are too large to treat by wave-function methods. The aim is not to review all that has been done using pathintegral dynamics methods, since there are already several excellent reviews in the literature [8][9][10][11], nor is it to cover allied topics, such as real-time path-integral methods [12][13][14][15][16][17] or the application of static path-integral methods to infer dynamical properties (e.g., tunnelling splittings [18,19] or quantum instanton rates [20]). Instead, I hope to convince a perhaps sceptical reader that (imaginary-time) path-integral dynamics methods are approximations to the exact quantum dynamics, which are necessarily crude when compared with wavefunction methods, but are often good enough to capture at least the essential physics. I will concentrate on single-surface Born-Oppenheimer (BO) dynamics, focusing mainly on the infrared spectrum of water.
To set the scene, let us write out the static pathintegral formalism for a one-dimensional system with hamiltonianĤ =p 2 /2m + V (q). The thermal expectation value 1 By 'quantum Boltzmann statistics' I mean the quantum statistics of distinguishable particles. a e-mail: sca10@cam.ac.uk (corresponding author) can be calculated using where A(q) N = 1 (2π ) N Z N dp dq A N (q) × e −β[HN (p,q)+SN (p,q)] , (with the quantum partition function Z N similarly defined) and where p ≡ {p 1 , . . . , p N }, q ≡ {q 1 , . . . , q N }, and β N = β/N . In other words, the trace over the diagonal matrix elements of the quantum Boltzmann operator e −βĤ is equivalent to a canonical phase-space average of an extended 'ring-polymer' system, consisting of N replicas or 'beads' of the original system, joined by harmonic springs. For a system such as liquid water at room temperature, N 30 is typically sufficient to give converged results [4] Efficient numerical methods have been devised to sample the Boltzmann ring-polymer distribution using Monte Carlo [2,5,21], or path-integral MD (PIMD) [3].
The early papers on PIMD warn the reader not to equate the ring-polymer trajectories to the real-time dynamics of the quantum system [3]. The exact quantum dynamics is of course given by the Schrödinger equation (or by real-time forward-backward Feynman paths) and the PIMD dynamics is used merely as a device for sampling the bead positions q l , with the bead momenta p l playing the role of dummy integration variables, with arbitrary masses (which have been set to the physical mass m in Eq. (4)). However, such was the desperation in the quantum dynamics community (where the application of wave-function methods in, say, 10 dimensions is even today a major challenge) that two heuristic methods were proposed for extracting realtime dynamical information from ring-polymer classical trajectories: centroid molecular dynamics (CMD) [23][24][25][26][27] and ring-polymer molecular dynamics (RPMD) [28][29][30][31][32].
In CMD, one propagates classical trajectories in the phase space of the centroid (ring-polymer centre of mass) on the potential of mean force obtained by integrating over the quantum thermal fluctuations (q l − Q 0 ). In RPMD, one propagates classical trajectories in the full 2N -dimensional phase space of the bead coordinates (p l , q l ), with the associated masses taken to be the physical masses [as they are in Eq. (4)]. One can show [23,28] that CMD and RPMD give the exact dynamics of the centroid in the high temperature (β → 0), zero time (t → 0), and harmonic limits, which ensures that they reproduce the exact quantum Kubo timecorrelation function (t.c.f.) of a pair of linear operators (i.e.,p orq) in these limits; RPMD is also exact at zero time for pairs of non-linear operators. A heuristic justification of CMD and RPMD is therefore that these methods interpolate the quantum Kubo t.c.f. between the limits just mentioned, and (by construction) they do so in a way that conserves the quantum Boltzmann distribution, which ensures that they correctly lock in zero-point energy and satisfy quantum detailed balance. However, there is more to these methods than mere interpolation between easy limits.
For example, CMD gives near quantitative agreement with experiment for the fundamental bands in the infrared spectrum of liquid water at room temperature [33,34]. Figure 1 shows (partially adiabatic) CMD predictions of these bands, computed using the accurate MB-pol potential energy [35,36] and dipole- Fig. 1 Experimental and simulated infrared spectrum of liquid water at 298.15 K. The simulations applied the partially adiabatic centroid molecular dynamics (PACMD) method to the MB-pol potential energy surface [35]. Reproduced with permission from Ref. [33]. Copyright 2015 American Chemical Society moment [33] surfaces. The agreement of the simulated fundamental bands with experiment is the closest yet achieved for a first-principles water simulation.
Another striking example is that RPMD is surprisingly accurate at predicting the rates of tunnelling reactions [29,30]. Viewed as an interpolation, one might expect RPMD to work for shallow tunnelling, where Wigner's correction factor [37] is valid, because the system penetrates only the parabolic tip of the reaction barrier; and also of course, the correct treatment by RPMD of zero-point energy is essential for quantum reaction rates [38,39]. However, it was found that RPMD works at far lower temperatures [29], where deep tunnelling speeds up the rate by orders of magnitude, even for the difficult case of highly asymmetric reaction barriers (for which an earlier proposed centroid transition-state theory [40] fails). In Ref. [41], it was shown that RPMD works so well for deep tunnelling, because the RPMD rate is controlled by a thermodynamic bottleneck consisting of thermal fluctuations of the ring-polymers around the instanton. 2 However, despite such successes, CMD and RPMD also fail to describe dynamics that one might expect to be not much more difficult than the scenarios just mentioned. Although CMD works well for water at room temperature, it suffers from a so-called 'curvature problem' as the temperature is lowered further [48,49], which artificially broadens and shifts to the red the OH-stretch band, making CMD inaccurate for ice [50]; RPMD avoids this problem but instead decorates the infrared spectra with spurious resonances [51,52] (caused by coupling of the polymer springs to the dynamics of the centroid), even at high temperatures.
This pattern of success and failure hinted that CMD and RPMD are approximations to an underlying theory that combines quantum (Boltzmann) statistics consistently with classical dynamics. That there should be such a 'classical dynamics-quantum statistics' theory was also consistent with earlier findings of more conventional semiclassical calculations [53], especially those based on the classical Wigner or 'linearized semiclassical-initial value representation' (LSC-IVR) method [53][54][55][56][57][58][59][60][61]. Such calculations gave strong evidence that most of the quantum effects in single-surface condensed-phase dynamics originate in the statistics. However, LSC-IVR can only be a short-time limit to such a theory, because it does not conserve the quantum Boltzmann distribution [52], and thus rapidly releases zero-point energy into inter-molecular degrees of freedom (e.g., causing ice to melt in less than a picosecond at 150 K) [62].
Reference [63] proposed a solution to this problem, which is that real-time coherence arises from 'jaggedness' in the imaginary-time Feynman paths. It is well known from simulations of static properties [5] that imaginary-time Feynman paths come in two varieties: jagged and smooth (Fig. 2). The jagged paths are the 'ring polymers' of Eqs. (3)-(7); such paths are never smooth, not even in the limit N → ∞. The smooth paths are obtained by Fourier-filtering the jagged paths. It is well known [5,64] that an ensemble of smooth paths converges to the same (exact) quantum thermal expectation values as an ensemble of jagged paths (albeit more slowly). However, Ref. [63] showed that this is not true for dynamics, as the effective Planck's constant associated with the smooth degrees of freedom is zero. Dynamics restricted to the smooth space is therefore classical. The statistics, however, is still fully quantum Boltzmann (because the smooth degrees freedom describe the quantum thermal fluctuations), and the dynamics, although classical, satisfies quantum detailed balance and thus correctly accounts for zero-point energy. This dynamics is called 'Matsubara dynamics' [63,[65][66][67][68][69][70].
These theoretical findings are backed up by numerical results, which show that, for simple models, Matsubara dynamics gives time-correlation functions which are in close agreement with the exact quantum results, indicating that the effects of real-time coherence are indeed small [66]. Moreover, it is easy to show that both CMD and RPMD are approximations to Matsubara dynamics, which result from taking respectively mean-field [66] and short-time limits [65]. However, a big disappointment with Matsubara dynamics is that it can- Fig. 2 Jagged 'ring-polymer' and smooth 'Matsubara' imaginary-time Feynman paths. The jaggedness is solely responsible for real-time quantum coherence; the dynamics of the smooth paths is therefore classical (and is known as 'Matsubara dynamics') [63] not be applied numerically to realistic systems, owing to a highly oscillatory phase factor [63]. Nonetheless, its relation with CMD and RPMD has turned out to be useful for explaining when these methods are valid, when and why they break down, and also for indicating how they can be improved (which has resulted in the recent QCMD method [71]).
The article is structured as follows. After reviewing jagged and smooth imaginary-time Feynman paths and their application to thermal expectation values (Sect. 2), we describe the chain of theoretical approximations that links the exact quantum dynamics of imaginary-time Feynman paths (Sect. 3) through Matsubara dynamics (Sect. 4) to mean-field Matsubara dynamics and CMD (Sect. 5) and (thermostatted) (T) RPMD (Sect. 6). We then review some of the numerical techniques used to make CMD and (T)RPMD practical (Sect. 7), describe the newly developed QCMD method (Sect. 8), and discuss some of the limitations of current path-integral methods (Sect. 9). Section 10 concludes the article.
Note Most equations will be written out for a onedimensional system, and the diagrams of Figs. 2, 3, and paths generated by the imaginary-time Schrödinger equation. It has long been known [5,64] that this jaggedness can be smoothed away using Fourier-filtering and that the resulting ensemble of smoothed paths gives the same static averages as the jagged paths. Starting with the 'jagged' ring-polymers of Eqs. (3)-(7), we can Fourier-filter them by converting to normal mode coordinates with N = (N − 1)/2 and 3 which diagonalise the spring term to where One can then smooth the paths by truncating |n| at some finite value M ≡ (M −1)/2, while taking N → ∞. Because |n| < M N , it follows from the form of T ln that the paths become smooth and differentiable functions of the imaginary time τ = 0 → β , which can be written: where the 'Matsubara frequencies' [72] are obtained by letting |n| N in Eq. (12). We will often refer to the modes Q n , |n| ≤ M as the 'Matsubara modes', the smooth space they inhabit as 'Matsubara space', and the smooth paths q(τ ) as 'Matsubara paths'. The lowest frequency mode Q 0 is the centroid of Eq. (8), and the modes Q n , n = 0 therefore describe the thermal quantum fluctuations around Q 0 (which vanish in the high temperature limit).
A similar smoothing of the (at present, fictitious) momentum variables p l to allows us to define the Matsubara thermal expectation value of a function A(q) of the position operator as where Z M is similarly defined, and with and P ≡ {P 0 , P ±1 , . . . , P ±M } and Q ≡ {Q 0 , Q ±1 , . . . , Q ±M }. Alternatively, we can write which regenerates Eqs. (17)-(18) on substituting for (p(τ ), q(τ )) using Eqs. (13)- (15). The function A M (Q) is defined similarly to U M (Q). If we make M large, the non-smooth M < |n| ≤ N modes correspond to very high frequencies ω n , and can therefore be integrated out analytically by steepestdescent [5,63], giving system-independent pre-factors which cancel in the numerator and denominator of Eq. (3). As a result i.e., the exact quantum thermal expectation value can be calculated using either smooth or jagged paths (although, in practice, the jagged paths give faster numerical convergence).

Momentum operators and phase-space path-integrals
It is straightforward to generalise the above to thermal expectation values of functions of the momentum operator. Inserting a Fourier identity, we can write which can then be evaluated using 'open' linear polymers [5,73] with the ends held apart at q ± Δ/2. To obtain the equivalent expression for smooth paths, it is better first to rewrite Eq. (23) more symmetrically as where which is identical to Eq. (23) since all p k , k = l can be integrated out. To smooth this expression, we convert to normal modes P n , Q n , take N → ∞ and truncate at |n| = M (as in Sect. 2.1), which gives where H M (P, Q) is given in Eq. (17), A M (P) is defined analogously to A M (Q), and As with A(q) , one can show that Equation (26) is a phase-space path-integral [74] which treatsp andq on an equal footing; unlike in Eq. (16), p(τ ) is no longer a dummy variable. Although the paths are smooth, they are also 'open', because there are no springs. Instead, there is a 'Matsubara phase' θ M (P, Q), which is what remains of the Fourier transforms in Eq. (24) after smoothing. The phase plays the equivalent role to the spring term S M (P, Q) of Eq. (16) in making the distribution quantum Boltzmann. To demonstrate this, we analytically continue [65] This expression shows that one can evaluate an expectation value over momentum operators using the identical smooth path-integral expression to that used for position operators in Eq. (16), provided that one uses A M (Π) in place of A M (P). It is easy to confirm that this expression is correct. For example, taking A(p) ≡p, we obtain p = lim M →∞ P 0 M = 0; taking A(p) to be the kinetic energy operatorT =p 2 /2m, we obtain which is the familiar 'primitive' kinetic energy estimator [22] expressed in terms of the Matsubara modes. Thus, smooth paths give the exact thermal expectation value for functions ofp as well as ofq. One can also show that the same is true for mixed functions ofp and q. However, major differences between the smooth and jagged spaces will appear when we consider dynamics in Sect. 4.

Exact ring-polymer quantum dynamics
Path-integral dynamics methods such as CMD and RPMD follow the time evolution of an entire imaginarytime Feynman path (or an ensemble of such paths); i.e., they follow the time evolution of a set of N replicas of the system. This way of describing dynamics seems at first sight to be at odds with the usual exact quantum description, which follows the time evolution of a single replica of the system. However, it has been known since a paper published by Shi and Geva [75] that one can in fact formulate exact quantum dynamics also as the dynamics of an entire imaginary-time Feynman path, and (as we will show in Sects. [4][5][6] this is the representation of the dynamics that is approximated by practical path-integral methods such as CMD and RPMD. 4 This is true provided AM (P) is an analytic function.
Consider first the more commonly used 'standard' version of the quantum t.c.f.

C [S]
which we can expand as where we assume for now thatÂ andB are functions of positionq (and we generalise to includep below). We can regard C AB (t) as a sum over Feynman loops; each loop has an 'open' imaginary-time section, connecting q ± Δ/2, which can be represented as a polymer 'string' (shown in blue in Fig. 3), and which is joined to realtime forward and backward sections (shown in red in Fig. 3). The latter represent the time evolution of a single replica of the system, which is located at q 1 at time zero and has evolved to z 1 by time t. The remaining replicas in the linear polymer (see Fig. 3) are therefore static and do not evolve in time; their job is to give the point q 1 a quantum Boltzmann weighting at time t = 0. Now, consider the Kubo-transformed version of the quantum t.c.f.
which contains equivalent physics to C in which G   AB (t). Following Shi and Geva [75], we insert N − 1 identities e iĤt/ e −iĤt/ ≡ 1 into the quantum Boltzmann operator, then let N → ∞, to obtain where and B N (z) is defined analogously to A N (q). This way of writing C AB (t) is equivalent to Eq. (35), but reveals two major differences between C showing that the integrand is invariant under cyclic permutations P (l → l + 1) of the bead index l. This property leads to the well-known [28] result that which is the same detailed balance relation satisfied by the classical t.c.f. This permutational invariance plays a key role in Matsubara dynamics (see Sect. 4). Second, the dynamics is no longer that of a single replica of the system, but instead that of N replicas, each located at one of the bead positions q l at t = 0, and evolving to z l at time t (see Fig. 3). In the limit N → ∞, each of the matrix elements q l−1 − Δ l−1 /2|e −βNĤ |q l + Δ l /2 can be replaced by a single ring-polymer spring, connecting q l−1 − Δ l−1 /2 to q l + Δ l /2. The dynamics of the replicas therefore corresponds to the exact quantum dynamics of a set of beads q l distributed at t=0 around a ring-polymer with N openings of width Δ l ; in other words, to the exact quantum dynamics of an imaginary-time Feynman path.
It is straightforward to generalise the form of C AB (t) given in Eq. (38) to momentum-dependentÂ andB. One inserts N Fourier identities into Eq. (35), and then follows the same steps that led to Eq. (38), obtaining [63] where and is the quantum Liouvillian [76]. Like Eq. (38), the integrand of Eq. (40) is invariant under P (l → l + 1), and the (exact quantum) dynamics propagated by e Lt is that of N replicas of the system, located initially at the points (p l , q l ) distributed around a phase-space Feynman path with N openings. Nobody would suggest using Eq. (40) to compute the quantum t.c.f., since it has increased the dimensionality of the system by a factor of N in the limit N → ∞. However, we show in the following sections that the exact quantum dynamics of the phase-space points (p l , q l ) is what lies behind pathintegral methods such as RPMD and CMD.

Matsubara dynamics
Matsubara dynamics emerges when the time-correlation function of Eq. (40) is smoothed. As mentioned in Sect. 2, smoothing has no effect on static properties (in the limit M → ∞). However, it has a major effect on the dynamics. Following similar steps to those of Sect. 2, the smoothed version of C AB (t) can be shown [63] to be where and H M (P, Q) and θ M (P, Q) are the Matsubara hamiltonian and phase defined in Eqs. (17)- (27). At t=0, C [Mats] AB (t) is exact; for the same reason that smoothed static expectation values are exact (Sect. 2). However, for t > 0, the dynamics in the smooth space is different from the dynamics in the jagged space. This is because, to avoid contaminating the dynamics with jagged modes (P n , Q n ), |n| > M , one needs to approximate the exact quantum Liouvillian L N of Eq. discarded [63]. When this is done, one notices that the non-classical terms in L M vanish, because they involve powers of N ≡ /N , in the limit N → ∞. As result, the Matsubara Liouvillian is simply and the dynamics is therefore classical in the space of the smooth modes (P n , Q n ). This is the 'Matsubara dynamics' referred to in the Introduction. We emphasise that the classicality has not been imposed a priori; it has arisen naturally, because the effective Planck's constant associated with the smooth space is zero. 5 All real-time coherence effects are therefore generated by the jagged modes (P n , Q n ), |n| > M . Smoothing out these modes to give Matsubara dynamics is thus equivalent to neglecting the effects of real-time coherence. 6 Since the dynamics in the smooth space is naturally classical, one would expect that it conserves the quantum Boltzmann distribution function e −βHM (P,Q) e −iβθM (P,Q) . This is indeed the case, and is easily shown to result from the invariance of the integrand of the (non-smoothed) C AB (t) under cyclic permutation P (l → l + 1) of the bead indices (see Sect. 2). After smoothing, this property survives as invariance under imaginary-time translation τ → τ + τ 0 , which is equivalent to an internal displacement of the smooth (p(τ ), q(τ )) path around itself (see Fig. 4). No torque can be exerted on this internal degree of freedom (since U M (Q) is invariant under τ → τ + τ 0 -see Eq. (19)), and therefore, the momentum conjugate to it is conserved. This momentum is (minus) the Matsubara phase θ M (P, Q) (as demonstrated in the Appendix). Writing out the multi-dimensional generalisation of the phase makes it clear that θ M (P, Q) is indeed a generalised internal angular momentum, with the components of p(τ ) tangential to q(τ ) integrated around the inside of the smooth loop. Since the hamiltonian H M (P, Q) is also conserved, it follows that Matsubara dynamics conserves the smoothed quantum Boltzmann distribution. Matsubara dynamics therefore gives a plausible explanation of how classical dynamics can emerge from a quantum Boltzmann distribution, without violating quantum detailed balance: the dynamics in such cases is dominated by the Matsubara modes (P n , Q n ), |n| ≤ M , which, being smooth, have an effective Planck's constant of zero. A growing body of evidence (see below) suggests that this picture is physically correct. However, Matsubara dynamics cannot be used as a practical simulation method without further approximation, because the phase θ M (P, Q) is usually too oscillatory to integrate over numerically. It is necessary to make further approximations to deal with θ M (P, Q), and these turn out to be the basis of CMD and RPMD [65].

Mean-field Matsubara dynamics and CMD
One can reduce the phase θ M (P, Q) by making M small, and replacing U M (Q) in Eq. (46) by the free energy obtained by integrating over the |n| > M modes (where N is a normalisation constant). This approach is referred to (perhaps confusingly) 7 as 'mean-field Matsubara dynamics' [66]. It is equivalent to approximating the exact quantum Liouvillian L N by its expectation value over the modes M < |n| ≤ N , in the limit N → ∞.
One could make mean-field Matsubara dynamics formally exact, using projection operators [77], but the  8 Reproduced from Ref. [71], with the permission of AIP Publishing resulting generalised Langevin equation would be difficult to interpret or implement (because of a complex friction kernel). Instead, we treat M as a convergence parameter. It seems reasonable to expect that the potential-of-mean-force dynamics will give a good approximation to the exact Matsubara dynamics when M is sufficiently large that the ring-polymer distribution around Q n , |n| ≤ M in Eq. (48) is compact.
Using this approach, it is possible to compute the linear-dipole infrared spectrum of a two-dimensional 'Champagne bottle' model of a rotating OH bond (with the t.c.f. multiplied by a damping function to mimic the decorrelation of liquid water) [66]. The results (shown in Fig. 5) are, to date, one of the few pieces of numerical evidence that directly confirm that Matsubara dynamics correctly accounts for the classical part of the dynamics in a quantum system at thermal equilibrium. 8 Note that the 200 K Matsubara spectrum dips below zero at about 260 cm −1 . This feature is an artifact of the truncation of the Fourier modes at M = 5 (see Ref. [73]).
The differences between the mean-field Matsubara and exact quantum results are small, namely a 20 cm −1 blue-shift and slight broadening of the stretch peak, which presumably result from the real-time quantum coherence associated with the neglected jagged modes (although a long convergence tail in M cannot be ruled out). Most importantly, the Matsubara OH-stretch frequency is independent of temperature, showing that it correctly accounts for zero-point energy.
Although the mean-field approach works for model systems, it does not reduce θ M (P, Q) sufficiently to treat realistic many-dimensional systems, unless one sets M = 1, which restricts the dynamics to the space of the centroid modes. In this case, the phase vanishes (since θ M (P, Q) is the momentum conjugate to motion around the centroid-see Fig. 4), and the dynamics reduces to classical dynamics on the free energy surface which is CMD [65,66]. As mentioned in the Introduction, CMD was first developed on the basis of heuristic arguments [23][24][25][26][27]. The above analysis tells us that, since CMD is meanfield Matsubara dynamics with M = 1, it is expected to work if the ring-polymer fluctuations around the centroid are compact, but to break down otherwise, when higher values of M will be required. For the twodimensional OH-stretch model, this compactness condition is satisfied perfectly at 800 K, less well at 600 and 400 K (c.f. the mean-field Matsubara and CMD results in Fig. 5), and breaks down badly at 200 K, where the OH-stretch band has a large artificial redshift and is broadened. This is an example of the socalled 'curvature problem' mentioned in the Introduction [48,49]. For the two-dimensional model, we were able to cure this problem by increasing M to give a compact ring-polymer distribution in Eq. (48) (which gave the mean-field Matsubara results of Fig. 5). However, for realistic systems, one cannot increase M > 1, and thus the CMD curvature problem can be a serious source of error at low temperatures (see Figs. 6 and 7). In Sect. 8, we explain how to use curvilinear coordinates to correct this problem using the recently developed QCMD method.

Ring-polymer molecular dynamics
Like CMD, RPMD was introduced on the basis of heuristic arguments [28][29][30][31][32], but can now be understood as an approximation to Matsubara dynamics [65]. Starting with the Matsubara t.c.f. of Eq. (45), we can analytically continue P n as in Eq. (30). This removes the phase, but analytically continues the Matsubara where and which propagates the dynamics into the complex plane, making the approach unusable. However, if we wish to compute only the short-time dynamics of the centroid, we can neglect the imaginary Liouvillian iL I , since it does not act directly on the centroid modes. The remaining real Liouvillian L RPMD describes smoothed RPMD, namely classical dynamics on the smoothed ring-polymer potential energy surface U M (Q)+S M (Q). Like static averages, the RPMD t.c.f. converges to the same results in the smooth space q(τ ) as in the jagged space q l , so practical RPMD calculations are always Fig. 7 The infrared spectrum of the q-TIP4P/F model of liquid water and ice, computed using QCMD, CMD, and TRPMD. Reproduced from Ref. [71], with the permission of AIP Publishing carried out in the jagged space (which yields faster numerical convergence).
As an approximation to Matsubara dynamics, RPMD is risky, because it is difficult to estimate how much time can elapse before the neglect of iL I starts to introduce noticeable errors into the dynamics of the centroid. Well-known symptoms of this failing are the spurious resonances which appear in RPMD simulations of infrared spectra [51,52], and which arise because the neglect of iL I causes the ring-polymer springs to shift the vibrational frequencies of the fluctuation modes, which in turn resonate with the dynamics of the centroid.
The other approach, which works for longer times, is thermostatted RPMD (TRPMD) [50,93,94]. Like plain RPMD, TRPMD was introduced on the basis of heuristic arguments [93]; it differs from RPMD only in attach-ing thermostats to the fluctuation modes Q n , n = 0 (leaving the centroid modes unthermostatted). The centroid therefore evolves unimpeded for as long as its dynamics is not affected by coupling to the fluctuation modes. When the coupling starts to become noticeable, the thermostat decorrelates the centroid dynamics, thus damping the t.c.f. As a result, TRPMD automatically follows the RPMD trajectory for as long as it gives a good approximation to the Matsubara dynamics of the centroid. This elegant ruse permits TRPMD to simulate infrared spectra, by eliminating the spurious resonances at the cost of artificially broadening the lineshapes. Figure 5 shows that, for the 2-dimensional OH model, TRPMD reproduces the position of the meanfield Matsubara OH-stretch frequencies (including the small ∼ 20 cm −1 blue-shift) across the entire temperature range, and that the artificial broadening grows with decreasing temperature (as the amplitudes of the fluctuation modes become larger); TRPMD behaves similarly for gas-phase and condensed-phase water (see Figs. 6 and 7 and Ref. [50]).

Practicalities of CMD and TRPMD simulations
In practical terms, CMD and TRPMD simulations are quite similar: classical trajectories are propagated on the ring-polymer potential energy surface V N (q) + S N (q), subject to a thermostat, which these days is usually a white-noise Langevin thermostat [97]. 10 For CMD simulations, this type of approach was developed as a method of computing the centroid force −∂F 0 (Q)/∂Q 0 on-the-fly [26,98,99]. The dynamics of the centroid are adiabatically separated from the dynamics of the n = 0 fluctuation modes, by making the masses of the latter much lighter than the (physical mass) of the centroid, and by applying a strong thermostat to the fluctuation modes (but not to the centroid). This adiabatic CMD (ACMD) approach makes it feasible to compute −∂F 0 (Q)/∂Q 0 on-the-fly, but the calculations can become costly, because the propagation time-step needs to be reduced to keep up with the rapidity of the fluctuation dynamics.
Accordingly, partially adiabatic (PA) CMD calculations are often performed to allow longer time-steps, in which the masses are reduced by less than is required to converge to the adiabatic limit [99]. The PACMD approach thus involves a certain amount of parameter tuning (of the mass-scaling and the thermostat strength), and can be thought of as an interpolation between RPMD (no mass-scaling, no thermostat) and CMD (adiabatic scaling, strong thermostat). In systems where the CMD curvature problem is mild, one can expect PACMD to give good results, in shifting the ring-polymer frequencies out of range of the spectrum (thus avoiding spurious resonances) [51], without picking up too much of a red shift and broaden- 10 Nosé-Hoover thermostats were used in earlier work.
ing. Such a regime is characterised by liquid water at 300 K [50,100,101]; PACMD calculations [33] were used to give the excellent agreement between MB-pol simulations and experiment shown in Fig. 1.
As mentioned above, TRPMD tends to give good spectral line positions (because it does not suffer from the curvature problem), but the thermostat (which is made sufficiently strong to damp off the spurious resonances) does significantly broaden the lineshapes [50,93]. If one can live with this broadening, TRPMD is definitely the path-integral method of choice for vibrational dynamics; it is cheaper than (PA)CMD, because the masses of the fluctuation modes are unscaled (permitting longer time steps), and it does not require parameter tuning (whereas in PACMD, both the adiabatic separation and the thermostat strength must be carefully chosen). Recent applications of TRPMD have included vibrational dynamics of aqueous proton defects [102], the solvated electron [103], and polyatomic organic molecules [104]. The TRPMD thermostat can also be coloured to reduce the broadening; this approach works well for the OH-stretch and bend bands in liquid water, but artificially blue-shifts the libration band [94].
Many of these TPRMD simulations (and also some PACMD and static simulations) used techniques based on ring-polymer contraction (RPC) [105][106][107][108][109] and multiple time-step (MTS) [109,110] propagation. The basis of these techniques is that the system potential can be decomposed into a cheap part, which needs to be evaluated at each polymer bead, and a costly part, which varies slowly around the ring-polymer and thus can be evaluated at fewer beads. A combination of RPC and MTS can speed up calculations by up to two orders of magnitude, thus allowing path-integral calculations to be combined with on-the-fly DFT potentials [102,109]. Another potential efficiency saving is the use of Cayley propagators [111,112], which also avoid possible nonergodicity artifacts. Closely related to TRPMD, the PIGLET approach [113] uses a strong quantum thermostat to reduce the number of polymer beads to just a few; it was previously used to compute static properties, but a post-processing approach has recently been developed which extracts dynamical properties from PIGLET, using an iterative image-space reconstruction algorithm [114]. The PIGLET vibrational density of states for liquid water is in reasonable agreement with the TRPMD result (see Fig. 8), although the stretch band is red shifted. Many of the above techniques are readily accessible to non-specialists via the i-PI pathintegral python interface [115].
We have focused mostly on vibrational spectra, but should point out that many of the above techniques can be used to improve the efficiency of RPMD calculations of thermal reaction rates [10,29,30,38,39,[82][83][84][85][86][87][88], and diffusion constants [31,32,[116][117][118][119][120][121][122] (for which which both CMD and TRPMD appear to work well, probably because any diffusive dynamics that satisfies quantum detailed balance and has the centroid as the diffusion coordinate is sufficient). Particularly noteworthy pathintegral simulations of diffusion constants include simu- The vibrational density of states of liquid water at 300 K, recovered from a PIGLET [113] simulation. Adapted from Ref. [114], with the permission of AIP Publishing lations of diffusion in quantum glasses [118], and in flexible q-TIP4P/F liquid water [117] (which revealed the importance of cancellation between competing intraand inter-molecular quantum effects). Recently, CMD has also been used to compute the thermal heat conductivity of liquid hydrogen and helium [123], using a reformulation of this observable that avoids non-linear Gibbs-Kubo relations [124].

Quasi-centroid molecular dynamics
As mentioned in Sect. 5, CMD gives a poor approximation to Matsubara dynamics at low temperatures, since the distribution of ring-polymers around the centroid becomes non-compact, spreading around the curve of the potential. When first identifying this problem, Marx and co-workers [48,49] suggested that it could be solved in principle by propagating the dynamics of curvilinear centroids, but the difficulty of formulating pathintegrals in curved spaces seems to have discouraged the development of such a method. Recently, however, Trenins et al. [71] have developed the quasicentroid molecular dynamics (QCMD) method, which uses curvilinear centroids to construct the potential of mean force, but propagates the dynamics in cartesians.
To illustrate QCMD, consider the two-dimensional OH-bond model discussed in Sect. 5. Figure 9 shows the ring-polymer distribution around the centroid during a CMD trajectory that suffers from the curvature problem. The distribution spreads around the curve of the potential, because the minimum of the centroid-constrained ring-polymer distribution has a delocalised 'string' geometry (shown in blue in Fig. 9) which describes an instantonic tunnelling path. Deep tunnelling is not expected to play a significant role in OHbond dynamics at 200 K, and instantonic paths clearly have a much lower Boltzmann weighting than compact paths at the same radial distance from the origin. The difficulty is that the centroid constraint has assigned these paths to the wrong radial distance, because the cartesian centroid lies at the 'focal point' of the curved instantonic paths, placing it at a smaller radial distance than any of the individual beads. As a result, the instantonic paths dominate and artificially lower the free energy at the inner turning point. To assign such paths to the correct radial distance (where they make a negligible contribution), one can constrain instead the polar centroids where r l = x 2 l + y 2 l and θ l = arctan(y l /x l ), and propagate the dynamics on the potential of mean force The ring-polymer distributions around the curvilinear centroids are then compact, as shown in Fig. 9. This compactness ensures that the dynamics on F 0 (R 0 ) gives a good approximation to the Matsubara dynamics of the polar centroid, and also that the polar centroid is close to the cartesian centroid (and is thus referred to as the 'quasi-centroid'), allowing it to be used in place of the cartesian centroid when computing time-correlation functions. This is the basis of the QCMD method. The compactness also ensures that rotation-vibration coupling between the fluctuation modes and the quasi-centroid is small, which in turn permits one to propagate the quasi-centroid trajectories in cartesian coordinates, invoking the polar coordinates only when generating F 0 (R 0 ). Figure 5 shows that the QCMD method yields excellent agreement with the mean-field Matsubara results for the vibrational spectrum of the two-dimensional OH-bond model. Using an adaptation of the adiabatic CMD approach (see Sect. 7), it has been possible recently to generalise QCMD to treat gas-phase water (which gives close agreement with the damped quantum results-see Fig. 6), and the q-TIP4P/F model of liquid water and ice (Fig. 7). The QCMD method thus has the poten- Fig. 9 Snapshots of CMD and QCMD ring-polymer bead distributions (red) taken from a simulation of the twodimensional OH-bond model at 200 K. At the inner turning point, the CMD distribution fluctuates around an instantonic geometry (thick blue curve), whereas the QCMD distribution fluctuates around a geometry at which all beads are collapsed to a point (blue dot). Reproduced from Ref. [71], with the permission of AIP Publishing tial to give accurate predictions of infrared spectra at lower temperatures than TRPMD and CMD. However, the QCMD method will need to be made more efficient before it can be applied to potential energy surfaces that are more expensive than the simple q-TIP4P/F model (as the adiabatic separation necessitates very small time-steps [71]), and the curvilinear centroid coordinates have to date not been generalised beyond water.

On the neglect of fluctuation dynamics
The CMD, TRPMD, and QCMD methods are alike in that they attempt to describe the dynamics of the centroid, but not that of the fluctuation modes [23,71,93]. So far, we have (tacitly) assumed that knowledge of the centroid dynamics is sufficient to determine the observables of interest. However, this is true only for time-correlation functions that involve pairs of linear operators. Non-linear operators also depend explicitly on the dynamics of the n = 0 fluctuation modes. For example, whenÂ =q 2 , A M (Q) = n Q 2 n . This does not necessarily mean that methods such as CMD and TRPMD break down for non-linear operators. However, it does mean that they can only be used to describe the centroid component of the operator. If the physical degrees of freedom in which the operator is non-linear involve, say, low-frequency librations, or long-range order parameters, then the centroid component is likely to be sufficient. This is probably the case for at least part of the non-linearity in the dipole-moment function of liquid water [33,125], which describes inhomogeneous variations in the local dipole moment between different hydrogen-bonded environments. However, for high-frequency, intramolecular, degrees of freedom, we can expect that the amplitudes of the fluctuation modes will be significant, such that their omission from the operator will introduce sizeable errrors.
It will be very difficult to overcome this limitation (except in the limit t → 0, where RPMD gives the exact quantum dynamics of any permutationally invariant function of the Q n [126]), since to propagate the fluctuation modes, one would need to sample the Matsubara phase. In Refs. [127][128][129], a planetary model was developed, which does capture the fluctuation dynamics to within a local harmonic approximation. This model gave good results for the dynamic structure factor of liquid hydrogen [128], and good overall predictions of the spectrum of liquid water [129], but it is too approximate to be used, for example, to investigate how nonlinearity in the dipole-moment surface of water affects the lineshapes.
A related problem is that the Matsubara dynamics of the fluctuation modes can also couple to the dynamics of the centroid. As mentioned in Sects. 5-8, practical path-integral methods neglect this coupling, under the assumption that the fluctuation dynamics can be either averaged out (CMD and QCMD) or damped (TRPMD). This assumption evidently works well for the fundamental bands of water (see Figs. 6, and 7), where the physical degrees of freedom are either highly anharmonic and classical (librations) or weakly anharmonic and quantum (vibrations). However, recent perturbation theory (PT) analyses [130,131] have shown that coupling of the fluctuation and centroid dynamics magnifies the amplitudes of the centroid overtone and combination vibrations by about an order of magnitude. Methods such as CMD, QCMD, and TRPMD therefore greatly underestimate overtone and combination band intensities. Again, the only way to correct rigorously for these deficiencies would be to include the dynamics of the fluctuation modes, which runs into the problem of the Matsubara phase. However, if one can assign the overtone and combination bands, one can at least correct the intensities approximately, using correction factors derived from PT [131]. The neglect of centroid-fluctuation coupling is also the most likely reason that CMD and TRPMD calculations of highly anharmonic vibrational frequencies give large blue-shifts [132] (although real-time coherence may also contribute to these errors).

Conclusions and recent developments
Path-integral dynamics methods have a rigorous basis in the exact quantum Kubo time-correlation function, which, upon smoothing of the imaginary-time Feynman paths, gives rise to a classical yet quantum-Boltzmannconserving dynamics called Matsubara dynamics [63]. Practical path-integral methods, such as CMD, (T)RPMD, and the recently developed QCMD, all make rather drastic approximations to Matsubara dynamics, but nevertheless manage to describe the dynamics of the centroid reasonably well, provided that it does not involve degrees of freedom which are simultaneously quantum and highly anharmonic. As a result, path-integral methods give a good approximation to, say, the fundamental bands in the spectrum of liquid water or ice, subject to the caveats mentioned above (namely, that CMD needs to be applied at sufficiently high temperatures to avoid the curvature problem, that TRPMD broadens the lineshapes, and that QCMD suffers from neither of these problems, but is currently 10-100 times more expensive than TRPMD).
In addition, RPMD gives an excellent approximation to the quantum rates of direct reactions (or indirect reactions that split into statistically decoupled direct steps), even when the rate is dominated by deep tunnelling, where it reproduces the instantaneous quantum flux through a bottleneck ensemble consisting of thermal fluctuations around the instanton. The CMD method also works in this way for the important special case of symmetric barriers. Both CMD and RPMD also give good results for diffusion coefficients.
Path-integral methods are expected to be costly, because they require the force to be evaluated along chains of replicas of the system. However, techniques such as RPC [105][106][107][108][109], MTS [109,110], and PIGLET [113,114] increasingly allow path-integral dynamics methods to be applied at not much more expense than a classical MD calculation, and the i-PI package [115] makes such techniques available to non-specialists.
As discussed throughout, path-integral calculations have some major limitations. It is very difficult to see how they can be generalised to treat non-linear operators when the non-linearity involves quantum degrees of freedom (except for the special case of instantaneous quantum fluxes in RPMD rate calculations), or to see how they can be made to accurately predict vibrational overtone intensities. For such properties, it might be better to use the LSC-IVR method [53][54][55][56][57][58][59][60][61], which, although failing to conserve the quantum Boltzmann distribution, does at least treat all the Matsubara modes on an equal footing (by following the dynamics of what is essentially a single polymer bead); another possible method might be the recently developed adaptive quantum thermal bath method of Plé et al. [130].
We finish by mentioning some areas of path-integral dynamics where progress has been made recently and which promises to extend the range of possible applications of path-integral methods.
Recent work has extended path-integral methods to non-equilibrium dynamics. The discussion of this article was restricted to dynamics at thermal equilibrium.
However, there appears to be no requirement when deriving Matsubara dynamics that the system be in thermal equilibrium: it is quite feasible to equilibrate, say, a system on one Born-Oppenheimer surface, before propagating it on another, or to equilibrate, say, parts of a system separately before combining them. All that is required is that the initial condition be describable in terms of smooth imaginary-time loops. This has been exploited recently with the development of nonequilibrium RPMD [133] (following earlier explorations of non-equilibrium CMD [134,135]), which has had interesting recent applications to hydrogen-graphene scattering [136,137], gas-surface reactions [138], and to excited-state dynamics [139]. Closely related to this are the recent developments of microcanonical RPMD [140][141][142].
Finally, much of the recent work on path-integral dynamics has looked at ways to extend these methods to treat non-adiabatic, multi-surface, dynamics [139,, building in part on earlier work [170][171][172]. The main challenge for non-adiabatic path-integral dynamics is the old problem of combining quantum (electronic) with classical (nuclear) degrees of freedom. Some calculations have used the surface-hopping technique to do this [139,146]. Another approach [147][148][149][150][151][152][153][154][155][156] is the use of mapping-variables [168,169] which has the advantage of making both the electronic and the nuclear degrees of freedom classical (at the price of requiring the use of projection operators). Recently, Chowdhury and Huo [70] have developed 'non-adiabatic Matsubara dynamics', which treats the mapping variables quantum mechanically and the nuclear coordinates by Matsubara dynamics. Progress has also been made in the treatment of charge-transfer reactions, with the development of Golden-rule ring-polymer transition-state approaches [157][158][159][160] and analytic continuations of instanton [161] and Wolynes theory [162,163,172], capable of treating the Marcus-inverted regime. However, it would be premature to say more at present about these rapidly developing areas: we leave this task to future review articles.
It is a pleasure to thank the group members, co-workers, and colleagues who have helped me to understand pathintegral dynamics over the last decade or so, many of whose works are referenced below. I owe additional thanks to George Trenins for reading through and commenting on the manuscript.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: All the numerical results discussed here have been published previously.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix: Equivalence of the Matsubara phase to a conserved internal angular momentum
To show that θM (P, Q) is equal to (minus) the momentum conjugate to τ0, we transform from the set of M Fourier modes Qn to a new set of coordinates, M − 1 of which (denoted ξ) are invariant under imaginary-time translation τ0, and one of which is τ0 itself [126]. The momentum conjugate to τ0 is then