We report the first numerical calculations in which converged Matsubara dynamics is compared directly with exact quantum dynamics with no artificial damping of the time-correlation functions (TCFs). The system treated is a Morse oscillator coupled to a harmonic bath. We show that, when the system–bath coupling is sufficiently strong, the Matsubara calculations can be converged by explicitly including up to *M* = 200 Matsubara modes, with the remaining modes included as a harmonic “tail” correction. The resulting Matsubara TCFs are in near-perfect agreement with the exact quantum TCFs, for non-linear as well as linear operators, at a temperature at which the TCFs are dominated by quantum thermal fluctuations. These results provide compelling evidence that incoherent classical dynamics can arise in the condensed phase at temperatures at which the statistics are dominated by quantum (Boltzmann) effects, as a result of smoothing of imaginary-time Feynman paths. The techniques developed here may also lead to efficient methods for benchmarking system–bath dynamics in the overdamped regime.

## I. INTRODUCTION

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^{34} and T-PIGS^{35} 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 $(\u223c100)$ 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 $(\u223c100)$ 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 $\u223c200$ 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.

## II. BACKGROUND THEORY

### A. Matsubara dynamics

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^{45}

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\u0304=(M\u22121)/2$ (so *M* is odd) $n\u0304\u2261\u2212n$, and *ω*_{n} = 2*πn*/*βℏ* are the Matsubara frequencies.^{46}

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\u0307t(\tau )=pt(\tau )/m$, $p\u0307t(\tau )=\u2212dV[qt(\tau )]/dqt(\tau )$. 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.^{47} The Matsubara dynamics approximation to the quantum Kubo TCF can then be shown^{1} 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 $A\u0302$ and $B\u0302$ in Eq. (1) are linear functions of $p\u0302$ and $q\u0302$, the functions $AP,Q$ and $BP(t),Q(t)$ in Eq. (4) become linear functions of the centroids *P*_{0} and *Q*_{0} [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^{39} 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 $A\u0302=B\u0302=q\u03022$, for which

### B. System–bath model

We consider the usual system–bath model^{48} 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,^{49}

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

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 $H\u0302$ and $H\u0302D$ 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 $H\u0302D$ means that Eq. (22) describes an artificial relaxation process and $C\u0303ABD(t)$ can be expected to be different to $C\u0303AB(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.

## III. SYSTEM–BATH MATSUBARA DYNAMICS

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).

### A. Matsubara GLE

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}(

**). Since we wish to approximate $C\u0303D(t)$ in Eq. (22), we also need the Matsubara Hamiltonian corresponding to**

*Q**H*

_{D}of Eq. (9), which is

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

where *∫*d*X*_{b} = *∏*_{α,n}*∫*d*X*_{αn} (and similarly for *P*_{b}) and

is the Matsubara phase. By analogy with the standard derivation of a classical GLE,^{48} 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).^{52} 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^{53}) by analytically continuing the bath modes, which amounts to replacing each *P*_{αn} by $P\alpha n+im\alpha \omega nX\alpha n\u0304$. Equation (25) then becomes

where

with $\omega \alpha n2=\omega \alpha 2+\omega 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.

### B. Analytically continuing the system variables

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\omega nQn\u0304$, which gives

where

and the equations of motion become^{54}

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*.

### C. Real-noise approximation

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

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^{56} for it to recover completely from the initial decorrelation of *X*_{αn} and $P\alpha n\u0304$, 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.

## IV. NUMERICAL TESTS OF THE MATSUBARA GLE

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 $\u27e8q\u03022q\u03022(t)\u27e9$.

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 $\u27e8q\u0302q\u0302(t)\u27e9$ and $\u27e8q\u03022q\u03022(t)\u27e9$, 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 $\u27e8q\u0302q\u0302(t)\u27e9$ 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.^{57}

System | D_{0} | D_{0,water}/2 |

D_{0,water} | 0.187 48 | |

ω_{well} | 1.703 04 × 10^{−2} | |

a | $\omega wellm/(2D0)$ | |

m | 1741.1 | |

Bath | η | 2η_{crit} = 2 × 2 mω_{well} |

ω_{c} | ω_{well} | |

Temperature | T | 150 K |

System | D_{0} | D_{0,water}/2 |

D_{0,water} | 0.187 48 | |

ω_{well} | 1.703 04 × 10^{−2} | |

a | $\omega wellm/(2D0)$ | |

m | 1741.1 | |

Bath | η | 2η_{crit} = 2 × 2 mω_{well} |

ω_{c} | ω_{well} | |

Temperature | T | 150 K |

However, the TCF we will focus on in what follows is the non-linear $\u27e8q\u03022q\u03022(t)\u27e9$. 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 $\u27e8q\u03022q\u03022(t)\u27e9$ 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 $\u27e8q\u03022q\u03022(t)\u27e9$ 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 $\u27e8q\u03022q\u03022(t)\u27e9$ 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 $\u27e8q\u03022q\u03022(t)\u27e9$.

### A. Matsubara GLE with complex noise

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^{5} trajectories being more than sufficient to converge the TCFs to graphical accuracy (and 2 × 10^{4} trajectories and just 100 implicit bath modes being sufficient to converge to within 1%). The resulting $\u27e8q\u03022q\u03022(t)\u27e9$ TCF for *M* = 25 is shown in Fig. 2. While the complex noise dynamics can be pushed up to *M* = 45, the resulting $\u27e8q\u03022q\u03022(t)\u27e9$ 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.

### B. Matsubara GLE with real noise

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*_{0} 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*.^{60} 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.

### C. Harmonic tail correction

Although *M* = 200 would appear to be a large number of Matsubara modes, the $\u27e8q\u03022q\u03022(t)\u27e9$ 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\u03022$ 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\u0304$. 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 $\u27e8q\u03022q\u03022(t)\u27e9$ 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].

### D. Agreement between Matsubara dynamics and HEOM

Figure 5 shows that the $\u27e8q\u03022q\u03022(t)\u27e9$ 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 $\u27e8q\u0302q\u0302(t)\u27e9$ 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.

## V. CONCLUSIONS

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 ($\u223c50$, 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.^{68}

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: POTENTIAL OF MEAN FORCE

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\u0304$ modes using an adiabatic decoupling parameter Γ as described in Ref. 39.

### APPENDIX B: EQUILIBRIUM MATSUBARA GLE

At thermal equilibrium, the Kubo-transformed TCF is

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

where $\zeta \u0302$ is the Laplace transform of the memory kernel

The analogous expression to Eq. (31) is then

where the driving terms are

and

which give Eqs. (35a) and (35b) of Sec. III A. When *∂K*_{n}(*t*)/*∂t* is defined, the two driving terms are related by^{69}

We then obtain

which is the equilibrium version of the Matsubara GLE [cf. the non-equilibrium direct-product GLE of Eq. (32)]. The noise term can be approximated to be real as described in Sec. III C. The driving terms *K*_{n} and *L*_{n} provide an additional memory, which means that even white noise simulations with a delta–function memory kernel are no longer Markovian.

### APPENDIX C: DEBYE BATH

The results given in Sec. IV used the Debye bath spectral density. Bath spectral densities are typically parameterized by the bath strength *η* (coupling strength) and cut-off frequency *ω*_{c} and are usually normalized such that

Our calculations used the following well known or easily derivable expressions for the Debye bath:

where we have used the efficient discretization of Ref. 59.

Corresponding expressions are easily derived for other spectral densities; see Ref. 70 for expressions for the white noise and Ohmic (with exponential cut-off) baths. In general, *K*_{n}, *L*_{n}, and $\zeta \u0302$ do not need to have closed analytic expressions, since they can easily be evaluated numerically.

### APPENDIX D: NUMERICAL INTEGRATION OF THE ANALYTICALLY CONTINUED MATSUBARA GLE

To integrate the Matsubara GLE, we employ the usual Liouvillian-splitting technique^{48,51} to obtain the propagator,

in which $m\omega V2q2/2$ is the harmonic part of *V*(*q*) and $\Pi n(t)=Pn(t)+im\omega nQn\u0304(t)$. Equation (D3) is a generalization of the “ring-polymer reference system propagator” of Ref. 71 to integrate the analytically continued equations of motion

### APPENDIX E: HARMONIC CORRECTION

To treat the modes $M\u0304<|n|\u2264Meff\u0304$ harmonically, we decompose the Matsubara TCF as

where $\u2211{Meff}$ denotes a sum over all the system Matsubara modes $|n|\u2264Meff\u0304$. We then further split the sums into

where *∑*^{{M}} denotes a sum over the $|n|\u2264M\u0304$ modes, which are included in the full anharmonic simulation, and *∑*^{{M}′^{}} a sum over the modes $M\u0304<|n|\u2264Meff\u0304$, which are treated harmonically. The purely harmonic terms are evaluated accurately from the numerical propagations carried out independently for each *n*, using $V(q)=m\omega V2q2/2$. The mixed terms are calculated using

which is exact if the modes *n*_{2} are harmonic. For sufficiently high |*n*| (typically |*n*| ≳ 200) the difference between the direct-product and equilibrium values was found to be negligible, allowing the $\u27e8Qn22(t)\u27e9$ terms to be replaced by their equilibrium values

Finally, Fig. 7 illustrates that the dynamics of the non-centroid Matsubara modes is largely harmonic.

## REFERENCES

Note that ω_{n} is negative for negative *n*.

Note that a different driving term of purely quantum origin is obtained for the equilibrium initial condition (see Appendix B).

Apart from making the substitution, one must move the limits of the integrals in Eq. (25) back onto the real line, which is justified by Cauchy’s residue theorem only if the integrand has no poles between the real axis and the line $P=\u2212im\omega nQn\u0304$. This condition is guaranteed at *t* = 0 and is likely at other times. A presence of a pole would lead to unstable trajectories, but the presence of unstable trajectories is not in itself a proof of the presence of a pole.

Presumably because the dynamics conserves the Matsubara phase.

This is the expected convergence behaviour, since the quantum Boltzmann weight can be expected to decrease monotonically for sufficiently large *M*. This trend is not obeyed once the trajectories become unstable, causing the integrals to diverge.

Note that we have not tested these approaches for anharmonic baths. We suspect that such systems would give rise to highly unstable Matsubara dynamics.

This is assuming that the derivative of *K*_{n}(*t*) is defined, which seems always to be the case except for a white noise bath at *t* = 0.