Comparison of Matsubara dynamics with exact quantum dynamics for an oscillator coupled to a dissipative bath

Matsubara dynamics was derived in Ref. 1 and put forward as a hypothesis to explain how classical dynamics can coexist with quantum Boltzmann statistics. Many approximate calculations over the years^2–21 have suggested that Born–Oppenheimer nuclear dynamics in the condensed phase often lives in or close to such a regime; for example, zero-point energy shifts in the infrared spectrum of liquid water are mainly corrected by including quantum (Boltzmann) effects into the statistics,^{17,18,22–24} and quantum rate coefficients can often be computed using just the quantum Boltzmann distribution at the barrier top.^{2,5,9,10,25–27} However, before Matsubara dynamics, no consistent theoretical description had been developed of such a “quantum Boltzmann–classical dynamics” regime. The nearest thing to it was the classical Wigner or LSC-IVR approach,^2–7 which is correct in the short-time limit, but fails at longer times because the plain Newtonian dynamics that it employs does not conserve the quantum Boltzmann distribution.

In Matsubara dynamics, by contrast, the classical dynamics takes place in an extended phase space obtained by describing the quantum Boltzmann statistics using imaginary-time Feynman paths. Crucially, these paths are smoothed. It has long been known^28–31 that smoothing the imaginary-time paths has no effect on static properties (although it does affect numerical convergence), but the unexpected finding in Ref. 1 was that the smoothing removes all the coherences from the real-time dynamics, because the effective Planck’s constant associated with the smooth space is zero. This property ensures that the Matsubara dynamics is classical but also quantum-Boltzmann conserving (without the need to add fictitious springs to the Hamiltonian). These properties suggest that Matsubara dynamics probably does correspond to the “quantum Boltzmann–classical dynamics” regime and useful progress has been made by making this assumption. In particular, the practical (thermostatted) ring-polymer molecular dynamics [(T)RPMD]^8–13 and centroid molecular dynamics (CMD)^15–18 methods, which were devised using heuristic arguments, have been shown to be (rather drastic) approximations to Matsubara dynamics; this understanding has led to the quasi-centroid molecular dynamics (QCMD) method,^19,32,33 which gives improved path-integral simulations of infrared spectra (and has led in turn to the powerful f-QCMD³⁴ and T-PIGS³⁵ methods). Matsubara dynamics has also been used to show that the well-known quantum enhancements of overtone intensities and Fermi splittings are caused mainly by the quantum Boltzmann statistics.^36,37

Given these and other potential uses of Matsubara dynamics, it would be good if one could compute it numerically for some model benchmark systems, in order to confirm that it does indeed correspond to the “quantum Boltzmann–classical dynamics” regime mentioned above. The difficulty is that Matsubara dynamics requires the imaginary-time paths to be represented as functions of position and momentum, which has the effect of converting the “polymer-springs” familiar from path-integral molecular dynamics^28,38 into a highly oscillatory phase factor. This “phase problem” has limited all previous tests of Matsubara dynamics to one- and two-dimensional models, to which it is not expected to apply, since there are far too few degrees of freedom for thermalization to decorrelate the time–correlation functions (TCFs) before real-time coherence effects become large. In Ref. 39, Trenins et al. sidestepped this problem for the simple “Champagne-bottle model” by artificially damping the TCFs. The resulting Matsubara-dynamics calculations of the infrared spectrum gave good but not perfect agreement with the (damped) exact quantum results, and the (small) errors were assumed to be caused by the real-time quantum coherence that survived the damping—but this could not be demonstrated conclusively. Another limitation of these calculations was that they could treat only TCFs involving linear operators, which depend explicitly only on the dynamics of the imaginary-time-path Cartesian centroids and thus require inclusion of only a few $(<10)$ non-centroid modes (describing the dynamical changes in shape of the paths). Many more $(∼100)$ non-centroid modes would be required to converge the Matsubara dynamics of non-linear operators, and one cannot take it for granted that Matsubara dynamics will be just as valid for non-linear operators as for linear ones, since all available approximate path-integral dynamics methods are known to fail for non-linear operators.

To confirm that Matsubara dynamics does correspond to the “quantum Boltzmann–classical dynamics” regime, one needs to apply it directly to a genuinely condensed-phase model in which the dynamics is expected to be mainly incoherent such that Matsubara dynamics should (if the hypothesis of Ref. 1 is correct) yield close agreement with the exact quantum results; and one needs to include large numbers $(∼100)$ of Matsubara modes to check that the dynamics is correct for non-linear operators. Recently, we had assumed that such a calculation would be computationally unrealistic, since the phase to be integrated over would be enormous. However, we show in this article that one can numerically converge the Matsubara dynamics for an anharmonic oscillator coupled to a harmonic bath, provided the system–bath coupling is sufficiently strong, in which case the resulting TCFs agree very closely with the exact quantum results, even for non-linear operators.

After summarizing Matsubara dynamics in Sec. II, we show in Sec. III how the Matsubara-dynamics approximation to an anharmonic oscillator coupled to a bath can be cast into the form of a generalized Langevin equation (GLE) with complex noise. In Sec. IV, we numerically show that the noise can be approximated to be real for a sufficiently strong system–bath coupling and that the real noise stabilizes the Matsubara dynamics when the system coordinates are analytically continued so as to convert the Matsubara phase into ring-polymer springs. Using this technique, we can compute the dynamics of $∼200$ Matsubara modes, which, with the use of a harmonic “tail” correction, is found to be sufficient to converge the Matsubara TCFs involving non-linear (position-squared) operators; these quantities agree very closely with the benchmark quantum results computed by solving hierarchical equations of motion (HEOM).^40–44 Section V concludes the article.

To derive Matsubara dynamics,^1,21,39 one starts with the exact quantum time-correlation function (TCF). The derivation follows most naturally from the Kubo-transformed TCF⁴⁵ 

where Z is the quantum partition function (and we will use Z in what follows as a general symbol to denote the partition function and any multiplicative constants that are required to normalize whatever distribution appears in the integral). One expresses the quantum dynamics as an (exact) propagation of the imaginary-time Feynman paths, then makes the approximation that these paths are smooth loops in imaginary time τ = 0 → βℏ. Keeping in one-dimension to simplify the notation (generalization to multi-dimensions is straightforward), these loops can be written as periodic functions p_t(τ), q_t(τ) of imaginary time. In practice, one expands q_t(τ) [and similarly p_t(τ)] as a Fourier series in terms of M Matsubara modes Q_n(t) as

where $M̄=(M−1)/2$ (so M is odd) $n̄≡−n$⁠, and ω_n = 2πn/βℏ are the Matsubara frequencies.⁴⁶ 

As first shown in Ref. 1, the effect of this smoothing is to strip out all real-time quantum coherences, so that the dynamics of (p_t(τ), q_t(τ)) becomes classical, with the equations of motion $q̇t(τ)=pt(τ)/m$⁠, $ṗt(τ)=−dV[qt(τ)]/dqt(τ)$⁠. In practice, one evaluates this dynamics in the space of the Matsubara modes (P_n(t), Q_n(t)) using the Hamiltonian

where the potential of mean force $UMQ$ is given in  Appendix A.⁴⁷ The Matsubara dynamics approximation to the quantum Kubo TCF can then be shown¹ to be the M → ∞ limit of

where

is the Matsubara phase.

The smoothness of (p_t(τ), q_t(τ)) ensures that the dynamics of (P_n(t), Q_n(t)) conserves the phase θ_M (because θ_M is the angular momentum conjugate to the internal rotation of the loop, on which there is no torque) and thus conserves the quantum Boltzmann distribution in Eq. (4). Unfortunately exp(−iβθ_M) is highly oscillatory, which makes it numerically extremely difficult to sample the distribution in Eq. (4); this is the “phase problem” mentioned in the Introduction.

When $Â$ and $B̂$ in Eq. (1) are linear functions of $p̂$ and $q̂$⁠, the functions $AP,Q$ and $BP(t),Q(t)$ in Eq. (4) become linear functions of the centroids P₀ and Q₀ [i.e., the centers of mass of the loops (p_t(τ), q_t(τ))]. This makes Matsubara dynamics much less numerically intractable (though still very difficult) for linear operators, since one only needs to know explicitly the Matsubara dynamics of the centroids, which often allows one to ignore the dynamics of all but a few of the n ≠ 0 modes describing the dynamical fluctuations around the centroids. This property was exploited by Trenins et al. when computing the spectrum of the champagne–bottle model³⁹ and is also the basis of the practical path-integral methods RPMD, CMD, and QCMD. An important goal of this article, however, is to deal with non-linear operators, which require that one calculates explicitly the dynamics of a large number of n ≠ 0 modes. In particular, we will treat $Â=B̂=q̂2$⁠, for which

We consider the usual system–bath model⁴⁸ with the classical Hamiltonian

where V(q) is the system potential and x_α, p_α and m_α are the positions, momenta, and masses of the bath oscillators, with frequency ω_α and coupling coefficients g_α. To facilitate a comparison with exact quantum benchmarks (see below), we calculate the direct-product TCF

in which the initial distribution is given by the direct-product Hamiltonian

but the Liouvillian $L$ is generated from H of Eq. (7) and therefore includes the system–bath coupling. The equations of motion for q(t) are solved in the same way as in the standard GLE derivation to obtain

which can be rewritten as the GLE

with the memory kernel

and random force term

where x_α, p_α, and q denote the values of these variables at initial time t = 0. The direct-product GLE of Eq. (11) is identical to the standard GLE except for the addition of the driving term −qζ(t), which acts as a short-term memory of the initial (direct-product) distribution.

It is customary to parameterize the bath via the spectral density

which determines both the memory kernel

and the random force, which are related by the fluctuation dissipation theorem,⁴⁹ 

Numerical solutions to Eq. (11) can be propagated using an explicit treatment of the bath modes or, equivalently, by discretizing the integrals in Eqs. (11) and (15) to give^50,51

where t_j ≡ jΔt, q_j ≡ q(t_j), and similarly for p_j and

where {w_α} are a set of quadrature weights. Comparing Eq. (18) with Eq. (12) gives

and

with

where λ_α and ξ_α are Gaussian random variates with a unit variance.

A variety of methods are available for treating the exact quantum version of the system–bath problem. Here, we will use the hierarchical equations of motion (HEOM) approach.^40–44 When applying HEOM, it is simpler (though not essential) to compute the direct-product Kubo TCF

in which $Ĥ$ and $ĤD$ are the quantum Hamiltonians corresponding to the full system–bath H of Eq. (7) and the direct-product system–bath H_D of Eq. (9), respectively. Clearly, the neglect of system–bath coupling in $ĤD$ means that Eq. (22) describes an artificial relaxation process and $C̃ABD(t)$ can be expected to be different to $C̃AB(t)$ (unless the bath is weak). This artificiality does not prevent us from using Eq. (22) to benchmark Matsubara dynamics; in fact, it is advantageous, since the thermalization depends critically on the relaxation dynamics of the non-centroid Matsubara modes (P_n, Q_n), making the comparison with Eq. (22) a stringent test.

In this section, we derive the Matsubara approximation to the direct-product Kubo TCF of Eq. (22) (Secs. III A and III B) and develop a real-noise approximation to it (Sec. III C). These derivations can easily be modified to give the Matsubara approximation to the equilibrium Kubo TCF of Eq. (1) (see  Appendix B).

The system–bath Matsubara Hamiltonian corresponding to H in Eq. (7) is easily shown to be

where the potential of mean force U_M(Q) is obtained by inserting the system potential V(q) into Eq. (A2) of  Appendix A. Note that each Matsubara mode Q_n is coupled to its own bath {X_αn} and {P_αn} and that the only term that couples Matsubara modes with different n is the system potential U_M(Q). Since we wish to approximate $C̃D(t)$ in Eq. (22), we also need the Matsubara Hamiltonian corresponding to H_D of Eq. (9), which is

The Matsubara approximation to the direct-product Kubo TCF of Eq. (22) can then be written

where ∫dX_b = ∏_α,n∫dX_αn (and similarly for P_b) and

is the Matsubara phase. By analogy with the standard derivation of a classical GLE,⁴⁸ it is straightforward to show that

where

The Matsubara equations of motion Eq. (27) resemble a generalization to M modes of the classical Eq. (11).⁵² However, G does not have the form of a random force term with Gaussian noise, because of the bath phase θ_B in Eq. (26). We, therefore, remove this phase (without making any approximation⁵³) by analytically continuing the bath modes, which amounts to replacing each P_αn by $Pαn+imαωnXαn̄$⁠. Equation (25) then becomes

where

with $ωαn2=ωα2+ωn2$ and the equations of motion retain the form of Eq. (27) except that G(Q, X_b, P_b) is replaced by

The variables P_b and X_b have Gaussian distributions in Eq. (29) and hence we can now write the equations of motion in the form of a GLE,

with a complex random force

where λ_αn and ξ_αn are Gaussian variates with unit variance. (Note that the random numbers in the imaginary parts of the sine coefficients are the same as those in the cosine coefficients of opposite n.)

We will refer to Eq. (32) as the Matsubara GLE. This equation looks superficially like a simple generalization to Matsubara modes of the classical GLE [Eq. (11)] and it reduces to CMD (classical dynamics on the centroid potential of mean force) when M = 1. However, when M > 1, the complex noise generated by Eq. (33) makes a major difference, since it satisfies a set of quantum fluctuation–dissipation relations

where

These relations ensure that the bath acts as a quantum thermostat such that the system Matsubara modes equilibrate to the (exact) quantum Boltzmann distribution. We emphasize that the only approximation made to obtain Eq. (32) from the exact quantum dynamics is to have smoothed the imaginary-time Feynman paths. Equation (32) is, therefore, exact, except for the neglect of real-time coherence.

The presence of the system Matsubara phase θ_S in Eq. (29) makes the integral very hard to integrate numerically (except when M = 1, which gives CMD, with θ_S = 0). We therefore convert θ_S to ring-polymer springs, by analytically continuing P_n, as was done for the bath modes in Sec. III A. This is equivalent to replacing each P_n in Eq. (29) by $Pn+imωnQn̄$⁠, which gives

where

and the equations of motion become⁵⁴ 

We emphasize that no additional approximation has been made here: Eqs. (36)–(38) are equivalent to Eqs. (25)–(28).

When the bathless versions of Eqs. (36)–(38) were derived in Ref. 54, it was pointed out that, although the analytic continuation has eliminated the phase, it has not eliminated the numerical difficulty because the dynamics now takes place in the complex plane and is found to be pathologically unstable (unless V is harmonic, M = 1, or the imaginary parts of the equations of motion are discarded to give RPMD). However, a surprising result reported in Sec. IV below is that the bath stabilizes the dynamics, up to some maximum value M, which depends on the bath strength and the anharmonicity of V.

One might expect that the stabilization just mentioned would be better if the random-force term were real since the imaginary part of this term will tend to kick the trajectories further along the imaginary axis. We, therefore, propose approximating the Matsubara GLE by

which is identical to Eq. (32), except that the imaginary part has been discarded from $RnC(t)$⁠, leaving the real part

Equations (36)–(38) remain the same, except that $RnC(t)$ is replaced by R_n(t) in Eq. (38).

To gauge how well this approximation is likely to work, we mention one point in its favor and one against. First, it was noted in Refs. 36 and 55 that one can eliminate the Matsubara phase by decorrelating the initial distribution in the TCF instead of by analytic continuation. When applied to the initial distribution of the bath modes

in Eq. (25), this “decorrelation approximation” gives

Then, if we retrace the steps in Secs. III A and III B (minus the analytic continuation of the bath modes), we obtain Eqs. (11) and (15). Hence, the neglect of the imaginary component of the noise is equivalent to approximating the initial distribution [by decorrelating the bath modes as in Eq. (42)], but it makes no approximation to the Matsubara dynamics equations of motion. We might, therefore, expect the errors due to neglecting $Im[RnC(t)]$ to be minor, since the neglect of correlations between the bath-modes will only take effect when bath modes with different n have had sufficient time to interact via the system potential V. However, the use of Eq. (42) also shows that the neglect of $Im[RnC(t)]$ modifies the quantum fluctuation–dissipation relations of Eq. (34) to

In other words, the Matsubara dynamics of the bath is insufficiently ergodic⁵⁶ for it to recover completely from the initial decorrelation of X_αn and $Pαn̄$⁠, which prevents it from correctly thermostatting to the exact quantum Boltzmann distribution. We can thus expect the real-noise approximation to the TCF not to tend to the exact t → ∞ limit. In Sec. IV, we find numerically that these errors decrease with increase in the bath strength and the stability gained by neglecting $Im[RnC(t)]$ is considerable.

In this Section, we test the stability of the analytically continued system Matsubara dynamics, for both complex and real noise, then investigate how far this takes us toward a numerically converged Matsubara-dynamics approximation to the direct-product Kubo TCF $⟨q̂2q̂2(t)⟩$⁠.

All calculations used a Morse system potential

and a Debye spectral density

for which expressions for bath properties such as g_αn, K_n, and L_n were derived (and are given in  Appendix C).

Figure 1 shows direct-product Kubo TCFs $⟨q̂q̂(t)⟩$ and $⟨q̂2q̂2(t)⟩$⁠, computed using the parameters of Table I, for a variety of standard methods. Evidently, this choice of parameters makes the dynamics significantly anharmonic. However, the anharmonicity causes only very weak coupling of the dynamics of the centroid to the dynamics of the n ≠ 0 Matsubara “fluctuation” modes, since the RPMD and CMD $⟨q̂q̂(t)⟩$ TCFs are in close agreement with HEOM, indicating that each of these methods gives a good description of the dynamics of the system centroid as it relaxes from the direct-product initial condition into equilibrium with the bath. These findings are consistent with previous applications of non-equilibrium RPMD to system–bath problems.⁵⁷ 

However, the TCF we will focus on in what follows is the non-linear $⟨q̂2q̂2(t)⟩$⁠. This TCF depends explicitly on the dynamics of the n ≠ 0 fluctuation modes, through Eq. (6), and it is well known^8,21 that RPMD and CMD do not describe this dynamics correctly. Figure 1 shows that the RPMD $⟨q̂2q̂2(t)⟩$ agrees with HEOM at t = 0 and at longer times, as expected (since it gives the correct quantum Boltzmann statistics in these limits), but is otherwise in poor agreement. The CMD results for $⟨q̂2q̂2(t)⟩$ are in extremely poor agreement with HEOM (and are not shown), for the obvious reason that CMD follows only the dynamics of the centroid, whereas most of the contribution to $⟨q̂2q̂2(t)⟩$ comes from the n ≠ 0 fluctuation modes. Our goal in what follows is to test numerically whether Matsubara dynamics is capable of describing the dynamics of these modes and thus giving good agreement with HEOM for $⟨q̂2q̂2(t)⟩$⁠.

As a first step, we tested the Matsubara GLE of Eq. (32) using the parameters of Table I. The system was equilibrated in a standard (one-dimensional) T(RPMD) simulation with N = 256 polymer beads, from which the analytically continued Matsubara GLE slave trajectories were propagated until t = 500 a.u. using the algorithm given in  Appendix D (with a timestep of 0.1 a.u. and including 1000 implicit bath modes). The Matsubara potential of mean force U_M was generated on the fly, following the procedure of Ref. 39, in which ring-polymer springs are attached to the N − M highest frequency ring-polymer normal modes, which follow adiabatically separated dynamics^16,58 subject to a strong PILE thermostat; an adiabaticity parameter of Γ = 16 (as defined in Ref. 39) was sufficient to converge the TCFs. The discretization of the bath frequencies was done using the efficient approach given in Ref. 59.

Surprisingly, we found that the coupling to the bath made the analytically continued system dynamics numerically stable for M ≤ 45. We also found that the number of trajectories needed to converge the TCFs was reasonable, with 2 × 10⁵ trajectories being more than sufficient to converge the TCFs to graphical accuracy (and 2 × 10⁴ trajectories and just 100 implicit bath modes being sufficient to converge to within 1%). The resulting $⟨q̂2q̂2(t)⟩$ TCF for M = 25 is shown in Fig. 2. While the complex noise dynamics can be pushed up to M = 45, the resulting $⟨q̂2q̂2(t)⟩$ TCF (not shown) still falls short of the exact quantum result shown in Fig. 1, indicating that many more Matsubara modes are required for numerical convergence.

We then repeated the calculation using the real-noise approximation given in Eqs. (39) and (40), keeping all other details of the calculation the same as in Sec. IV A. As discussed in Sec. III C, the real-noise approximation is not guaranteed to thermalize correctly to the exact quantum Boltzmann distribution, because it erroneously replaces the fluctuation–dissipation relations of Eq. (34) by Eq. (43). However, Fig. 2 shows that the resulting thermostatting error is small and that it decreases with an increase in η such that η = 2η_crit gives errors that are within graphical accuracy for M = 25 and the error for η = η_crit is still very small (within 1.6%). Thus, the imaginary component of the noise only has a subtle effect on the equilibration dynamics once η exceeds a certain value (which, in this system, is about η = η_crit); future studies will be needed to determine why this is so.

Very usefully, the real-noise approximation was found to make the analytically continued Matsubara dynamics much more stable, increasing the maximum value of M from M = 45 (complex noise) to M ≈ 200 (real noise). Since the stability is likely to depend on the bath strength and the anharmonicity of V(q), we varied the bath parameters η and ω_c and the system parameter D₀ to investigate the maximum value of M at which less than 0.1% of the trajectories diverged. This criterion was used because it was found to ensure that the gaps between TCFs computed using M and M + 2 modes decreased on increasing M.⁶⁰ Figure 3 shows that, as expected, stability increases with an increasing bath strength (i.e., increasing either η or ω_c) and decreases with increasing anharmonicity in V(q). Note that, even for η < η_crit, a considerable number of Matsubara modes can still be propagated stably.

Although M = 200 would appear to be a large number of Matsubara modes, the $⟨q̂2q̂2(t)⟩$ TCF requires a much larger value of M for numerical convergence. This is on account of the slowness with which the explicit sum over the Matsubara modes in Eq. (6) converges with respect to M, as has been noted in the context of static path integral calculations by Doll and co-workers.^28–30 We can gauge the likely length of this “Matsubara tail” by computing the (static) thermal expectation value of $q̂2$ at 150 K for a (bathless) harmonic oscillator with frequency ω_well. Figure 4 shows that thousands of Matsubara modes are needed to converge to a result attainable with just 100 ring-polymer beads. This is the price one pays for smoothing the imaginary-time Feynman paths.

However, Doll and co-workers showed^28–30 that it is possible to converge smooth path-integrals to the exact static quantities by treating all the Matsubara modes harmonically above some large value of |n| (which is treated as a convergence parameter). We therefore expect that the Matsubara dynamics of these modes can also be treated harmonically, since the highly oscillatory Matsubara phases will confine the dynamics of Q_n to a small region around zero. This is borne out by Fig. 7 of  Appendix E, which shows that the centroid dynamics is much more anharmonic than that of the n ≠ 0 fluctuation modes, which quickly become harmonic with increasing n.^61,72,73 In the harmonic limit, the Matsubara GLE splits into pairs of GLEs, each of which couples only the modes n and $n̄$⁠. These equations are straightforward to solve numerically, since the analytically continued harmonic trajectories are stable.

We, therefore, solved the analytically continued (real-noise) Matsubara GLE [Eq. (39)] for M modes, then generated harmonic corrections for the |n| > (M − 1)/2 modes up to a large value |n| = (M_eff − 1)/2 (as described in  Appendix E). For the system given in Table I, the tail correction converged the $⟨q̂2q̂2(t)⟩$ TCF using M = 55 and M_eff = 10 001, as demonstrated in Fig. 5 [where identical curves are obtained within graphical accuracy for M = 55 and M = 85, demonstrating that the dynamics of the |n| > (M − 1)/2 Matsubara modes is harmonic].

Figure 5 shows that the $⟨q̂2q̂2(t)⟩$ TCF computed using the real-noise Matsubara GLE for η = 2η_crit is in almost perfect agreement with the numerically exact HEOM results. We also found that the linear $⟨q̂q̂(t)⟩$ real-noise Matsubara TCF (not shown) agrees with the HEOM results within graphical accuracy and has converged by M = 15 (thus, correcting the very small discrepancy between CMD and HEOM in Fig. 1).

Figure 6 shows what happens when we apply the Matsubara GLE at weaker bath-strengths. For η = η_crit the Matsubara TCF is in very close agreement with the HEOM result, with a barely noticeable drift at longer times. For η = 0.5η_crit, the Matsubara and HEOM TCFs are also in close agreement until about t = 250 a.u., after which the errors from the real-noise approximation (cf. Fig. 2) cause a significant drift in the thermalization.

Finally, we also plot in Fig. 6 (see also Fig. 5(b) the result of truncating the calculation at M = 1, with all modes 0 < |n| < (M_eff − 1)/2 included using the harmonic correction of Sec. IV C. This approach is equivalent to using CMD to describe the motion of the centroid and treating the fluctuation dynamics harmonically. These results still have significant errors, but they are a marked improvement on the purely harmonic result (see Fig. 1 for η = 2η_crit); this is consistent with the findings given in Sec. IV C that the dynamics of the centroid is much more anharmonic than that of the fluctuations.

The results given in Sec. IV are the first converged Matsubara dynamics TCFs to be compared with exact quantum TCFs (without artificially damping in time) and to be obtained for non-linear operators, which report on the dynamics of the entire delocalized quantum Boltzmann distribution. The almost perfect agreement between the Matsubara dynamics and the exact quantum results for the two highest bath strengths confirms that Matsubara dynamics corresponds to the “quantum Boltzmann–classical dynamics” regime, and that this regime arises because thermal averaging decouples the dynamics of the smooth imaginary-time Feynman paths from the jagged paths.

The use of Matsubara dynamics also gives fresh physical insight into this old (system–bath) problem. It has long been known that the system–bath coupling enhances the effective strength of the “ring-polymer springs” along the system coordinate, causing delocalized static features, such as instantons, to become more compact.^62,63 The results given in Sec. IV show that the dynamics behaves analogously, with the system fluctuation modes (which correspond to the “spring” degrees of freedom) behaving much more harmonically than the centroid. Nonetheless, sufficient anharmonicity remains that a considerable number of fluctuation modes (⁠$∼50$⁠, in the calculations of Sec. IV) need to be coupled to the dynamics of the centroid to give quantitative agreement with the exact quantum result.

This near-harmonic dynamics of the fluctuation modes is probably what made the Matsubara calculations numerically tractable. We found that we could analytically continue the dynamics into the complex plane (thus converting the Matsubara phase into springs), where it remained numerically stable up to some maximum number of modes M, which increased with the bath coupling strength and decreased with the system anharmonicity. This suggests that Matsubara dynamics may be useful as a practical method for treating system–bath problems in the overdamped limit, where exact quantum methods such as HEOM^40–44 and QUAPI^64–66 (which are very efficient for treating relatively weak system–bath dampings) are known to become expensive. So far, we have treated dynamics on a single Born–Oppenheimer surface, but perhaps approaches such as that given in Ref. 67 can be used to extend the treatment developed here to non-adiabatic dynamics.⁶⁸ 

The authors thank Jimin Li, Robert Jack, Ilija Srpak, and William Moore for extensive discussions about the work in this paper. A.P. acknowledges an EPSRC DTA Ph.D. studentship from the UK Engineering and Physical Sciences Research Council; E.S.P. acknowledges an External Research Studentship from Trinity College, Cambridge.

The authors have no conflicts to disclose.

Adam Prada: Conceptualization (supporting); Formal analysis (equal); Investigation (lead); Methodology (lead); Software (lead); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). Eszter S. Pós: Conceptualization (supporting); Investigation (supporting); Methodology (supporting); Writing – review & editing (supporting). Stuart C. Althorpe: Conceptualization (lead); Formal analysis (equal); Investigation (supporting); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Following Ref. 39, the potential of mean force is defined by

with

and

We evaluated U_M(Q) on the fly by artificially increasing the frequencies of the $|n|>M̄$ modes using an adiabatic decoupling parameter Γ as described in Ref. 39.

At thermal equilibrium, the Kubo-transformed TCF is

Following the procedure given in Sec. III A, we obtain

where $ζ̂$ is the Laplace transform of the memory kernel