An exact representation of quantum mechanics using the language of phase-space variables provides a natural starting point to introduce and develop semiclassical approximations for the calculation of time correlation functions. Here, we introduce an exact path-integral formalism for calculations of multi-time quantum correlation functions as canonical averages over ring-polymer dynamics in imaginary time. The formulation provides a general formalism that exploits the symmetry of path integrals with respect to permutations in imaginary time, expressing correlations as products of imaginary-time-translation-invariant phase-space functions coupled through Poisson bracket operators. The method naturally recovers the classical limit of multi-time correlation functions and provides an interpretation of quantum dynamics in terms of “interfering trajectories” of the ring-polymer in phase space. The introduced phase-space formulation provides a rigorous framework for the future development of quantum dynamics methods that exploit the invariance of imaginary time path integrals to cyclic permutations.

Quantum thermal time correlation functions (TCFs) are ubiquitous in chemistry and physics and play a central role in the description of the dynamical properties of condensed phase systems.1–3 Indeed, (multi-)TCFs are the key constituent elements of dynamical theories involving linear and nonlinear response,4,5 chemical dynamics,2 and even quantum chaos.6 However, the exact full quantum mechanical calculations of TCFs for condensed phase systems involving hundreds of degrees of freedom are still impractical. As such, the development of approximate methods that capture nuclear quantum effects—such as tunneling and zero point energy—while relying on (semi)classical dynamics is of great interest to the chemistry and physics community.

The representation of quantum mechanics based on the language of phase-space variables provides a natural starting point for the implementation of semiclassical approximations. For equilibrium properties, the Feynman path-integral representation of quantum mechanics7,8 provides an appealing exact framework, allowing for calculations based on classical phase-space dynamics in imaginary time. In that framework, the canonical Boltzmann distribution of a quantum particle becomes isomorphic to the phase-space distribution of a classical ring-polymer composed of particle replica “beads” linked by harmonic springs, as defined by the discretized imaginary-time Feynman paths.9 Thermal quantum averages are thus obtained as canonical phase-space integrals over the ring-polymer ensemble. The ensemble of ring-polymer configurations is generated using classical methods, such as Monte-Carlo sampling or molecular dynamics simulations, thus allowing for simulations of systems in the condensed phase.9–12 

Extension of imaginary time path-integral techniques to calculate dynamical properties is more challenging due to the presence of the time evolution propagator controlling the dynamical evolution. Although one can express the dynamical evolution in terms of real-time forward–backward Feynman paths,13 the resulting time correlation functions involve phase factors that are difficult to evaluate numerically. (See Refs. 14–16 for recent developments to address that challenge.) As such, a variety of approximate methods that follow the time evolution of an entire imaginary time path, such as Matsubara dynamics,17 centroid molecular dynamics (CMD),18,19 and ring-polymer molecular dynamics (RPMD),20 have been proposed and applied with relative success for the calculation of one-time correlation functions. Understanding how these methods can be derived from approximations of the exact quantum Liouvillian and, especially, how to apply and extend these techniques to calculations of multi-time correlation functions is a subject of great research interest.21–23,28

Here, we introduce a formally exact formulation of multi-time correlation functions in terms of imaginary time Feynman paths with dynamics controlled by the exact ring-polymer quantum Liouvillian. The formulation provides a general formalism that exploits the symmetry of path integrals with respect to permutations in imaginary time, expressing correlations as products of imaginary-time-translation-invariant phase-space functions coupled through Poisson bracket operators and providing an interpretation of quantum dynamics in terms of “interfering trajectories” of the ring-polymer in phase space. Moreover, the framework reduces to classical multi-time correlation functions when → 0, providing a natural derivation of the classical limit of quantum mechanics.

The resulting phase-space formulation establishes a rigorous and appealing framework for the development of novel theoretical tools for studying quantum dynamics effects in condensed phase systems. Indeed, since the exact quantum Liouvillian governs the time evolution of ring-polymer observables, the full quantum-mechanical evaluation of the quantum correlation functions is impractical but for the simplest system. The introduction of semiclassical approximations is, therefore, needed to effectively evaluate the multi-time correlation functions. Remarkably, the invariance of the ring-polymer to imaginary time translation has been recently used to derive quantum Boltzmann-preserving semiclassical approximations to single-time17,18,24–26 and two-times22,27 Kubo transformed TCFs. Similar approximations can be applied to the generalized TCFs presented in this work, providing a variety of path-integral-based semiclassical approximations for the evaluation of multi-time correlation functions. (A detailed derivation of these approximations will be provided elsewhere.) The introduced phase-space formulation serves as the starting point for the development of novel semiclassical approximations.

The paper is organized as follows: Sec. II introduces the notation for calculations of equilibrium quantum averages as imaginary time path integrals in phase space. Section III extends the ring-polymer phase-space representation for calculations of multi-time correlation functions. Section IV discusses the properties of the imaginary time path-integrals as well as the symmetries and properties of ring-polymer multi-time correlation functions. Section V summarizes and concludes.

In this section, we introduce the notation and key elements for the imaginary-time path-integral phase-space representation of quantum mechanics. We remark that several elements of this section were originally introduced in Refs. 17 and 21. To keep the notation simple, we will focus on a one-dimensional system with Hamiltonian Ĥ=p̂2/2m+V(q̂). Nevertheless, the resulting expressions are also applicable to higher-dimensional systems.

We focus on the thermal ensemble average of an operator Ô=O(q̂,p̂),

Ô=Z1TreβĤÔ,
(1)

where

Z=TreβĤ
(2)

is the partition function and β = 1/kBT, with T the temperature and kB the Boltzmann constant.

We introduce the ring-polymer phase-space average,17 

ÔNZN1(2π)NdqdpeβĤN̄(q,p)ÔN(q,p),
(3)

to compute expectation values in terms of phase-space averages with ∫dq∫dq1∫dqN. Here, the N ring-polymer coordinates q = {q1, …, qN} are defined by the discretized imaginary time Feynman path subject to cyclic boundary conditions (i.e., q0 = qN). Analogous definitions apply for the ring-polymer of momenta p = {p1, …, pN}.

The ring-polymer phase-space function ÔN(q,p), introduced by Eq. (3), is the generalized Wigner–Weyl transform defined, as follows:43 

ÔN(q,p)1Nj=1NÔW(qj,pj),
(4)

where

ÔW(q,p)dΔeipΔqΔ2Ôq+Δ2,
(5)

is the Wigner–Weyl transform of the operator Ô,29–32 providing a map between operators in Hilbert space and classical-like functions in phase space for the ring-polymer variables (q, p) (see Fig. 1). For simplicity, we will often suppress the (q, p) dependence. We remark that for operators that depend only on position,

O(q̂)N(q,p)=1Nj=1NO(qj),
(6)

so the generalized Wigner–Weyl transform reduces to a ring-polymer average.9 Similar considerations hold for operators that only depend on momentum.

FIG. 1.

(Left) Schematic representation of the generalized Boltzmann–Wigner transform [eβĤ]N̄(q,p) [Eq. (7)]. Blue segments represent imaginary time propagation y|eβNĤ|x. (Right) Schematic representation of the generalized Wigner transform [Ô]N(q,p) [Eq. (4)]. The initial ring-polymer distribution [Ô]N(q,p), associated with observable Ô, consists of an average over single-bead terms.

FIG. 1.

(Left) Schematic representation of the generalized Boltzmann–Wigner transform [eβĤ]N̄(q,p) [Eq. (7)]. Blue segments represent imaginary time propagation y|eβNĤ|x. (Right) Schematic representation of the generalized Wigner transform [Ô]N(q,p) [Eq. (4)]. The initial ring-polymer distribution [Ô]N(q,p), associated with observable Ô, consists of an average over single-bead terms.

Close modal

The (pseudo)-probability distribution of the ring-polymer phase-space configuration (q, p), introduced by Eq. (3), is given by the generalized Boltzmann–Wigner transform eβĤN̄(q,p) defined, as follows:17,21,28

eβĤN̄(q,p)dΔl=1NeiplΔlql1Δl12eβNĤql+Δl2,
(7)

with βN = β/N and normalization constant given by

ZN1(2π)NdqdpeβĤN̄(q,p).
(8)

By evaluating the matrix elements y|eβNĤ|x with the symmetric Trotter approximation in the N → ∞ limit, an explicit form for eβĤN̄ can be obtained (see Sec. I C of the supplementary material) as follows:

eβĤN̄=m2πβN2N/2dΔeS(q,p,Δ),
(9)

with

S(q,p,Δ)=il=1NplΔlm2βN2l=1Nql+Δl2ql1+Δl122βN2l=1NVql+Δl2+VqlΔl2.
(10)

Note that Eq. (9) represents a ring-polymer comprised of N replicas of the system corresponding to the “beads” with coordinates ql distributed around an open ring-polymer with N openings of fixed width Δl33 (see Fig. 1). Note that upon integration over momenta p, the generalized Boltzmann–Wigner transform reduces to

dpeβĤN̄=2πmβNN/2eβNU(q),
(11)

with

U(q)=m2βN22l=1N(qlql1)2+l=1NV(ql),
(12)

which corresponds to the standard path-integral expression involving a close ring-polymer of beads connected with harmonic springs between nearest-neighbors.9 

The significance of the phase-space average introduced by Eq. (3) is that it provides an exact representation of quantum mechanical ensemble averages. In fact, evaluating the integrals over the Wigner transform, we obtain (see Secs. I A and I B of the supplementary material)

Ô=ÔN.
(13)

Therefore, Eq. (3) does not involve any kind of approximation, being equivalent to Eq. (1). It is just a generalization of the standard Wigner–Weyl phase-space representation29–32 to ring-polymer variables, allowing for numerically exact evaluations of quantum ensemble averages as canonical phase-space averages of an extended ring-polymer system.7–9 

Dynamical extensions of Eqs. (3) and (4) can be obtained for time-dependent operators, Ô(t)=eiĤt/ÔeiĤt/. In fact, the ring-polymer phase-space representation introduced in Sec. II A remains valid when introducing the substitution ÔÔ(t) in Eq. (4) giving the time-dependent phase-space distribution [Ô(t)]N(q,p). We remark that the Hamiltonian involved in the time evolution operator e±iĤt/ could differ from the one governing the statistics eβĤ, allowing the formalism introduced in this work to be applied to systems out of equilibrium.

The phase-space time evolution of [Ô(t)]N(q,p) is obtained by extending the standard Liouvillian formalism34,35 to ring-polymer coordinates, a formulation that is also convenient for the development of semiclassical approximations.

We introduce the Janus ring-polymer operator,21 

ΛNj=1NpjqjqjpjΛj,
(14)

where the arrows indicate the direction in which the derivative is applied. Note that ΛN is defined as minus the Poisson bracket operator for the phase-space variables of the ring-polymer. We also introduce the sine and cosine coupling operators s and c, which are defined as follows:

s2Nsin2ΛN,
(15)

and

ccos2ΛN.
(16)

 Appendix A summarizes some useful properties of s and c. Finally, we introduce the ring-polymer quantum Liouvillian operator,17,21,28

LNĤNs=j=1Npjmqj2V(qj)sin2qjpjLj,
(17)

which represents a generalization of the Moyal expansion of the quantum Liouvillian29,36 to ring-polymer phase-space variables. We remark that the potential energy V in Eq. (17) describes the dynamics of the system (not necessarily its statistics). Note that the quantum Liouvillian can be expressed in powers of 2 by Taylor expansion of the sine function, and that for harmonic potentials, for which higher-order derivatives vanish, the classical Liouvillian is recovered.

The Liouvillian formulation allows for a compact notation and interpretation of time-evolved observables in ring-polymer phase space (see Sec. I E in supplementary material), namely21,24

Ô(t)N(q,p)=eLNtÔN(q,p)=j=1NeLjtÔW(qj,pj).
(18)

In this sense, one can evaluate time-dependent averages, as follows:

Ô(t)N=ZN1(2π)NdqdpeβĤN̄(q,p)eLNtÔN(q,p).
(19)

For systems at equilibrium, for which eLNteβĤN̄=eβĤN̄, Eq. (19) reduces to the time-independent ensemble average, introduced by Eq. (3). We note that Eq. (19) is formally exact, allowing for numerically exact quantum dynamics simulations based on LN as averages over the contributions of ring-polymer “trajectories” in phase space. In effect, the initial ring-polymer phase-space distribution [Ô]N(q,p), which involves the sum over single-beads terms [Eq. (4)], is evolved by the propagator eLNt, according to the exact quantum Liouvillian LN, describing the motion of N independent replicas of the system17,33,37 (Fig. 2). When the potential energy is harmonic, the resulting propagation is analogous to classical propagation, albeit in an extended phase space. However, when the potential is anharmonic, the quantum Liouvillian includes higher-order derivatives beyond the classical terms so,

eLNt[Ô]N(q,p)[Ô]N(eLNtq,eLNtp),
(20)

and the amplitude of the Wigner transform at phase-space point (q(t), p(t)) [with q(t) and p(t) obtained by classical evolution of q(0) and p(0)] no longer corresponds to the time-evolved Wigner transform [Ô(t)]N(q,p). Only for potentials with vanishing higher-order derivatives (i.e., harmonic potential or free particle), Eq. (20) becomes an equality and the ring-polymer follows a classical trajectory in phase space. Nevertheless, as shown in Sec. III, thinking in terms of quantum trajectories provides a useful interpretation of multipoint time correlation functions.

FIG. 2.

Schematic representation of the time evolution of ring-polymer phase-space distribution [Ô(t)]N(q,p) [Eq. (18)]. The initial ring-polymer distribution [Ô]N(q,p) associated with observable Ô, which consists of single-bead terms, is evolved in time by the exact quantum Liouvillian LN [Eq. (17)], which evolves individual replicas of the system.

FIG. 2.

Schematic representation of the time evolution of ring-polymer phase-space distribution [Ô(t)]N(q,p) [Eq. (18)]. The initial ring-polymer distribution [Ô]N(q,p) associated with observable Ô, which consists of single-bead terms, is evolved in time by the exact quantum Liouvillian LN [Eq. (17)], which evolves individual replicas of the system.

Close modal

Section II introduced a path-integral phase-space representation of quantum mechanics that allows for the calculation of ensemble or one-time-dependent averages. However, many dynamical properties of condensed phase systems are encoded in multi-time correlation functions. In this section, we combine the ring-polymer statistics [eβĤ]N̄ with the time evolution described by the ring-polymer propagator eLNt[Ô]N to derive and introduce an exact ring-polymer phase-space representation of multipoint time correlation functions.

We introduce the “sine” and “cosine” two-point correlation functions in ring-polymer phase space as follows:

B̂(t1)sÂ(t0)NZN1(2π)NdqdpeβĤN̄×B̂(t1)NsÂ(t0)N
(21)

and

B̂(t1)cÂ(t0)NZN1(2π)NdqdpeβĤN̄×B̂(t1)NcÂ(t0)N,
(22)

where [Ô(t)]N=eLNt[Ô]N represents a quantum ring-polymer trajectory [Eq. (18)] and we understand that in the previous and subsequent equations, the sine/cosine operators couple pairs of time-evolved observables. The resulting sine and cosine correlation functions given by Eqs. (21) and (22) have an intuitive interpretation in terms of phase-space trajectories of the ring-polymer (Fig. 3): The two initial ring-polymer phase-space functions [Â]N and [B̂]N are evolved in time by the propagator eLNt and coupled by a sine or cosine couplings operator.

FIG. 3.

Schematic representation of the ring-polymer phase-space two-point time correlation functions. The initial ring-polymer distributions [Ô]N(q,p) associated with observables  and B̂ [Eq. (4)] evolve in time by the exact quantum Liouvillian LN (see Fig. 2). The sine/cosine coupling operators s/c [Eqs. (15) and (16)] couple observables at different times, giving rise to different quantum correlation functions [Eqs. (21) and (22)].

FIG. 3.

Schematic representation of the ring-polymer phase-space two-point time correlation functions. The initial ring-polymer distributions [Ô]N(q,p) associated with observables  and B̂ [Eq. (4)] evolve in time by the exact quantum Liouvillian LN (see Fig. 2). The sine/cosine coupling operators s/c [Eqs. (15) and (16)] couple observables at different times, giving rise to different quantum correlation functions [Eqs. (21) and (22)].

Close modal

The significance of the sine/cosine ring-polymer correlation functions introduced by Eqs. (21) and (22) is that they have a one-to-one correspondence to correlation functions defined in Hilbert space. Indeed, in  Appendix B, we prove that the sine correlation function satisfies

B̂(t1)sÂ(t0)N=iB̂(t1),Â(t0),
(23)

with [Ô1,Ô0]=Ô1Ô0Ô0Ô1 denoting a commutator. On the other hand, in  Appendix B, we show that the cosine correlation function is related to the Kubo transformed (KT) time correlation function given by

limNB̂(t1)cÂ(t0)N=B̂(t1)Â(t0),
(24)

where we define the KT correlation function of two arbitrary operators as38 

Ô1Ô2=1β0βdλÔ1(iλ)Ô2.
(25)

(See  Appendix C for properties of the KT.) In other words, in the limit N → ∞ Eqs. (21) and (22) give an exact path-integral phase-space representation of time correlation functions that, according to Eqs. (23) and (24), involve a commutation relation or a Kubo integral between two observables Â(t0) and B̂(t1). This is the first main result of this paper.44 

The relevance of the sine and cosine two-point correlation functions is clear since the commutator correlation function [Eq. (23)] and the Kubo transform correlation function [Eq. (24)] are related through the quantum fluctuation–dissipation relation,1,2

iB̂(t1),Â(t0)=βddt0B̂(t1)Â(t0),
(26)

and are central to linear response theory ( Appendix D). Therefore, the sine and cosine functions allow for ring-polymer calculations of transport coefficients, reaction rates, and simulations of (linear) spectroscopy.5,38

The path-integral phase-space formulation can be extended to calculations of multi-time correlation functions. Here, we consider the “sine–sine” three-point correlation function, defined in ring-polymer phase space, as follows:

Ĉ(t2)sB̂(t1)sÂ(t0)NZN1(2π)NdqdpeβĤN̄×Ĉ(t2)NsB̂(t1)NsÂ(t0)N.
(27)

Figure 4 shows a schematic representation of this correlation corresponding to the evolution of three observables Â, B̂, and Ĉ coupled by sine coupling operators. Note that since in our notation sine and cosine operators couple pairs of functions, for the case of correlation functions involving more than two operators one needs to specify in which order the operators are coupled [see Eq. (A3)].

FIG. 4.

Schematic representation of the ring-polymer phase-space three-point time correlation functions [Eqs. (27), (29) and (30)]. The ring-polymer distributions [Ô]N(q,p) associated with observables Â, B̂, and Ĉ [Eq. (4)] evolve in time by the exact quantum Liouvillian LN (Fig. 2). The sine/cosine coupling operators s/c [Eqs. (15) and (16)] couple observables at different times, giving rise to different quantum correlation functions [Eqs. (27), (29) and (30)].

FIG. 4.

Schematic representation of the ring-polymer phase-space three-point time correlation functions [Eqs. (27), (29) and (30)]. The ring-polymer distributions [Ô]N(q,p) associated with observables Â, B̂, and Ĉ [Eq. (4)] evolve in time by the exact quantum Liouvillian LN (Fig. 2). The sine/cosine coupling operators s/c [Eqs. (15) and (16)] couple observables at different times, giving rise to different quantum correlation functions [Eqs. (27), (29) and (30)].

Close modal

The significance of Eq. (27) is that it has a correspondence with a three-point correlation in Hilbert space. Specifically, in  Appendix B, we show that

Ĉ(t2)sB̂(t1)sÂ(t0)N=i2Ĉ(t2),B̂(t1),Â(t0),
(28)

so, the “sine–sine” correlation function is an exact path-integral representation of a correlation function involving a double commutation relation between observables Â(t0), B̂(t1), and Ĉ(t2). To the best of our knowledge, Eq. (28) provides a novel representation of a multi-time correlation function and represents the first main result of this section.

Remarkably, similar considerations apply to the “sine–cosine” and “cosine–cosine” correlation functions,

Ĉ(t2)sB̂(t1)cÂ(t0)NZN1(2π)NdqdpeβĤN̄×Ĉ(t2)NsB̂(t1)NcÂ(t0)N,
(29)

and

Ĉ(t2)cB̂(t1)cÂ(t0)NZN1(2π)NdqdpeβĤN̄×Ĉ(t2)NcB̂(t1)NcÂ(t0)N.
(30)

(We remark that for ring-polymer correlation functions involving only the cosine coupling operator, the order in which the operators are coupled is irrelevant [see Eq. (A3c)].) Specifically, in  Appendix B, we prove that the “cosine–cosine” correlation function provides a ring-polymer phase-space representation to the symmetrized Double Kubo transformed (DKT) time-correlation function,

limNĈ(t2)cB̂(t1)cÂ(t0)N=Ĉ(t2)B̂(t1)Â(t0),
(31)

where we define the symmetrized DKT time correlation function as21,22,27,39

Ô1Ô2Ô3=1β20βdλ00βdλ1T̂βÔ1(iλ0)Ô2(iλ1)Ô3,
(32)

with T̂β an imaginary time-ordering operator that orders the product of operators so their imaginary time arguments increase from right to left and ensure that there is no backward imaginary time propagation inside the integral [Eq. (C5)].21 (See  Appendix C for properties of the DKT.) On the other hand, the “sine–cosine” correlation function provides a phase-space representation of a KT quantum correlation function involving a commutator relation, namely,

limNĈ(t2)sB̂(t1)cÂ(t0)N=iĈ(t2),B̂(t1)Â(t0).
(33)

The ring-polymer phase-space correlation functions defined by Eqs. (27) and (29) and (30), and their connection to the exact quantum correlation functions [Eqs. (28), (31) and (33)] represent the second main result of this paper.

It is interesting that all the Hilbert space correlation functions defined by Eqs. (28), (31) and (33) naturally appear in the context of second-order response theory (see  Appendix D). As such, the three-point sine/cosine correlation functions are relevant for ring-polymer calculations of second-order spectroscopy.

The previous sections and  Appendix B showed that the path-integral phase-space representation can be generalized to calculations of two-point and three-point time correlation functions. Generalizations to multipoint time correlation functions are also possible, although the number of combinations of sine/cosine couplings between time-evolved observables quickly grows with the order of the correlation and the derivation of the expressions become tedious. Here, we introduce a convenient mapping between the expressions of correlations in ring-polymer phase space and correlations in Hilbert space.

The relation introduced by Eq. (23) suggests that the sine coupling operator is related to the commutator as follows:

Ô1NsÔ2NiÔ1,Ô2,
(34)

while Eq. (24) connects the cosine coupling operator with the Kubo integral, as follows:

Ô1NcÔ2N1β0βdλÔ1(iλ)Ô2=Ô1Ô2.
(35)

With the mapping introduced by Eqs. (34) and (35), the connection between canonical averages in ring-polymer phase space and traces in Hilbert space [Eqs. (23) and (24)] immediately follow. More interestingly, these associations still hold for the case of three-point correlation functions [see Eqs. (28), (31) and (33)] as long as one respects the ordering of operations of the sine coupling [e.g., Eq. (28)] and adopt the convention that correlation functions involving multiple Kubo integrals should be symmetrized in imaginary time [e.g., Eq. (31)]. With those considerations, Eqs. (34) and (35) provide a mapping between operators in Hilbert space and their corresponding ring-polymer representations, allowing to extend the path-integral phase-space formulation for general multipoint correlation functions.

A general n-point ring-polymer phase-space correlation function can be defined as follows:

Ôn(tn)JnJ2Ô1(t1)J1Ô0(t0)NZN1(2π)NdqdpeβĤN̄×Ôn(tn)NJnJ2Ô1(t1)NJ1Ô0(t0)N,
(36)

where each Jj, j = 1, …, n, represents either a sine coupling operator (s) or a cosine coupling operator (c). Note that since sine and cosine operators couple pairs of functions, one needs to specify a particular associative rule in the previous equation by defining which pairs of observables are coupled to each other and in which order (see below for examples). Note that different multipoint time correlation functions can be generated by changing the order of the correlation (n), the type of coupling between observables (s/c), and the particular associative rule. For n = 1 and n = 2, one recovers the correlations defined in Secs. III A and III B, respectively. We remark that Eq. (36) represents the most general exact ring-polymer phase-space representation of a time correlation function and, when combined with the associations introduced by Eqs. (34) and (35) mapping it to Hilbert space correlations, represents the main result of this paper.

As an example of how to apply the present theory, let us consider the n = 3 case and focus on the correlation involving three sine coupling operators. Using the mapping rule of Eq. (34), we postulate the equality

D̂(t3)sĈ(t2)sB̂(t1)sÂ(t0)N=i3D̂(t3),Ĉ(t2),B̂(t1),Â(t0)
(37)

between ring-polymer and Hilbert space correlation functions. By extending the derivation in  Appendix B, it is not difficult to show that the previous equality holds, proving that the ring-polymer correlation defined by the left-hand side of the expression above is an exact path-integral phase-space representation of the correlation defined by the right-hand side.

As a second example, consider the four-point correlations involving two sine and one cosine couplings. Using the mapping rules of Eqs. (34) and (35), we postulate the equivalences

D̂(t3)sĈ(t2)sB̂(t1)cÂ(t0)N=i2D̂(t3),Ĉ(t2),B̂(t1)Â(t0)
(38)

and

D̂(t3)sĈ(t2)cB̂(t1)sÂ(t0)N=i2D̂(t3),Ĉ(t2)B̂(t1),Â(t0)
(39)

between ring-polymer functions and KT correlation functions. Note that the correlations presented above exemplify two different association rules between observables giving rise to different Hilbert space functions (i.e., a KT involving one double commutator vs two single commutators). A straightforward extension of the derivation in  Appendix B shows that the previous equalities hold in the N → ∞ limit, validating the mapping rules Eqs. (34) and (35). Similar considerations follow for other types of couplings (although the math becomes cumbersome).

Incidentally, note that the four-point correlations functions defined above naturally appear in the context of third-order response theory (see  Appendix D). As such, the ring-polymer phase-space representation of multipoint time correlation functions provides an exact path-integral formulation of response theory.

The ring-polymer phase-space representation of general multipoint time correlation functions provided by Eq. (36) [and for two- and three-point functions, in particular, by Eqs. (21), (22), (27), (29) and (30)] represent the main result of this paper. The present formulation is exact and, in combination with the mapping rules introduced by Eqs. (34) and (35), provides a route to evaluate Hilbert space correlation functions in terms of ring-polymer phase-space averages that, to the best of our knowledge, represent a novel result. In this section, we discuss several interesting properties of the phase-space formulation.

The ring-polymer phase-space correlation functions are symmetric with respect to imaginary time translation, i.e., they are invariant to cyclic permutations of coordinates of the path integral.17,24 In fact, note that both the Boltzmann factor [eβĤ]N̄ and Wigner–Weyl functions [Ô]N, as well as the Liouvillian LN and Janus operator ΛN, are invariant to the cyclic permutation qlql+1. As such, the ring-polymer phase space provides a representation that emphasizes this symmetry for different multipoint correlation functions. Note that the invariance to translation in imaginary time is a well-known symmetry of the standard Kubo transformed correlation functions17,21 and has been recently used for the derivation of quantum Boltzmann-preserving semiclassical approximations for the calculation of KT and DKT correlation functions.17,21,22,24,25,28 To the best of our knowledge, this is the first time that a general expression that emphasizes this symmetry for general multipoint TCF has been presented.

The ring-polymer phase-space representation highlights the symmetry and the similarity between the different time correlation functions in ring-polymer phase space (Figs. 3 and 4). In particular, notice that each commutation relation in Hilbert space maps into a sine coupling operator (s) in ring-polymer phase space, whereas the presence of a Kubo integral contributes with a cosine coupling operator (c). In this sense, the different time correlation functions in ring-polymer space only differ in the type of coupling between observables. Note that if not for the sine and cosine coupling terms, the time evolution of the correlation functions would involve the action of the exact Liouvillian operator LN on the classical-like observables [O]N at different times independently (see Figs. 3 and 4). The effect of the sine/cosine terms is to couple and correlate the observables at different times with one another, incorporating effects due to interference and noncommutative in the correlation function, giving rise to different quantum correlation functions.

The structure of the ring-polymer phase-space functions [Eq. (36)] also provides an interesting interpretation of the quantum correlations in terms of “interfering trajectories”40 albeit in the extended ring-polymer phase space. In effect, note that since

sn=02n(ΛN)2n+1,
(40)
cn=02n(ΛN)2n,
(41)

the ring-polymer phase-space functions can be expressed in terms of an infinite series in powers of of correlations involving terms of the form [O1(t1)]N(ΛN)n[O0(t0)]N (with allowed values of n according to the sine/cosine parity). Using the definition of the ring-polymer Janus operator [Eq. (14)], each term in the expansion can be expressed as a hierarchy of different order stability matrices of the form

Mijk(n)(t)=n[O(t)]Nxi(0)xj(0)xk(0),
(42)

with xi = (qi, pi) denoting a ring-polymer phase-space point. Note that in classical mechanics, the stability matrix describes the sensitivity of a trajectory at time t to changes in the initial conditions and naturally appears in the context of classical response theory and classical chaos.40–42 Equations (36) and (42) naturally generalized this concept to the ring-polymer phase space. We remark that since the ring-polymer observables are propagated according to the exact Liouvillian operator and the hierarchy of stability matrices Mijk(n)(t) contribute to infinite orders in in the sine/cosine expansions, Eq. (36) is quantum-mechanically exact. This interpretation of the ring-polymer formulation provides a rigorous framework for the development of semiclassical approximations.

Due to the properties of the sine and cosine operators, there are different but equivalent forms of expressing time correlation functions in ring-polymer phase space. Depending on the problem at hand, some notations could be more advantageous than others, allowing to simplify the problem. For example, one could use relations Eqs. (A2a) and (A2b) to express the correlations in terms of [eβĤ]N̄c[Â(t0)]N or [eβĤ]N̄s[Â(t0)]N factors instead. Note that if [Â(t0)]N is only first order in position or momenta, only the first term in the cosine/sine expansion will survive and [eβĤ]N̄c[Â(t0)]N = [eβĤ]N̄[Â(t0)]N, whereas [eβĤ]N̄s[Â(t0)]N[eβĤ]N̄ΛN[Â(t0)]N. We remark that previously derived expressions for the KT and DKT by others17,28 and us21 are expressed in these alternative forms (see Sec. III of the supplementary material). The notation used in this paper, however, provides the most general form of correlation functions, allowing to apply the theory to general (nonlinear) operators and nonequilibrium systems and highlighting the similarities and symmetries of the different correlation functions.

Before concluding, it is instructive to analyze the emergence of the classical limit from the ring-polymer phase-space representation. Classical analogs of the general correlation functions defined by Eq. (36) can be obtained by reducing the phase space to one bead (i.e., setting N = 1) and noting that

lim0c=1,
(43a)
lim0s=Λ,
(43b)
lim0L=pmqV(q)qp,
(43c)

along with the fact that the quantum Boltzmann distribution [Eq. (7)] can be replaced with the classical distribution in the → 0 limit. For example, in the case of two-point correlations, the sine and cosine TCFs reduce to

B̂(t1)sÂ(t0)N0N=1B(t1)ΛA(t0)cl,
(44a)
B̂(t1)cÂ(t0)N0N=1B(t1)A(t0)cl,
(44b)

whereas cl represents an average over the classical Boltzmann distribution eβH, the observables O(t) are now classical dynamical variables that evolve following classical trajectories, and O1ΛO2{O1,O2}PB represents the negative of the classical Poisson bracket. In the case of three-point correlations, one obtains

Ĉ(t2)sB̂(t1)sÂ(t0)N0N=1C(t2)ΛB(t1)ΛA(t0)cl,
(45a)
Ĉ(t2)sB̂(t1)cÂ(t0)N0N=1C(t2)ΛB(t1)A(t0)cl,
(45b)
Ĉ(t2)cB̂(t1)cÂ(t0)N0N=1C(t2)B(t1)A(t0)cl.
(45c)

The classical expressions presented in Eqs. (44) and (45), and the more general expressions that can be derived from Eq. (36) in the classical limit, correspond to TCF that naturally appears in the context of classical response theory and classical chaos theory.40–42 As such, the exact ring-polymer phase-space representation presented in Sec. III establishes a rigorous and attractive framework that provides a bridge between the quantum and classical limits, allowing for the development of novel theoretical/computational tools for the incorporation of nuclear quantum effects in multi-time correlation functions.

In this work, we have presented a general and exact quantum dynamical formulation of multi-time correlation functions in terms of classical-like canonical averages over an extended ring-polymer phase space [Eq. (36)]. The ring-polymer representation provides a general formalism that exploits the symmetry of path integrals with respect to permutations in imaginary time, expressing correlations as products of imaginary-time-translation-invariant phase-space functions coupled through Poisson bracket operators. Moreover, the formulation provides a natural interpretation of quantum dynamics in terms of “interfering trajectories,” albeit in an extended phase space, and recovers the classical limit of multi-time correlation functions when → 0.

The ring-polymer phase-space formulation of multipoint time correlation functions presented in this work is formally exact. However, the full quantum-mechanical evaluation of the quantum correlation functions is impractical but for the simplest systems. For example, the time evolution of ring-polymer phase-space distributions [Eq. (18)] involves the exact Liouvillian and is tantamount to solving the Schrödinger equation. On the other hand, the sine/cosine coupling operators in this formulation involve the evaluation of ΛN interactions to (in principle) infinite order, although for operators that can be expressed as polynomials in ring-polymer phase-space coordinates, the expansion is exactly truncated at some order. The application and development of semiclassical approximations are, therefore, needed to effectively evaluate multi-time correlation functions. The present ring-polymer formulation provides a rigorous and attractive framework for the development of novel theoretical/computational tools for studying quantum dynamics effects in large/complex molecular systems.

The supplementary material includes details of the derivation of the imaginary time path-integral phase-space formulation.

P.E.V. thanks Kenneth A. Jung and Daniel Laria for insightful discussions. V.S.B. acknowledges support from the NSF Grant No. CHE-1900160 and high-performance computing time from the Yale High Performance Computing Center.

The authors have no conflicts to disclose.

Pablo E. Videla: Conceptualization (equal); Investigation (lead); Writing – original draft (equal); Writing – review & editing (equal). Victor S. Batista: Conceptualization (equal); Funding acquisition (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

It is worth highlighting some properties of the sine and cosine coupling operators (see Sec. I G of the supplementary material). The sine and cosine coupling operators are anti-self-adjoint and self-adjoint, respectively, i.e.,

Ô1NsÔ2N=Ô2NsÔ1N,
(A1a)
Ô1NcÔ2N=Ô2NcÔ1N.
(A1b)

Moreover, considering the integral over the extended ring-polymer phase space, it is straightforward to verify that the following properties hold:

dqdpeβĤN̄Ô1NsÔ2N=dqdpeβĤN̄sÔ1NÔ2N,
(A2a)
dqdpeβĤN̄Ô1NcÔ2N=dqdpeβĤN̄cÔ1NÔ2N.
(A2b)

The sine and cosine coupling operators applied to multiple generalized Wigner transforms satisfied the following associative rules:

Ô1NsÔ2NsÔ3N=Ô1NsÔ2NsÔ3N+Ô1NsÔ3NsÔ2N,
(A3a)
Ô1NcÔ2NsÔ3N=Ô1NcÔ2NsÔ3N+Ô1NsÔ3NcÔ2N,
(A3b)
Ô1NcÔ2NcÔ3N=Ô1NcÔ2NcÔ3N.
(A3c)

Note that Eq. (A3a) represents the Jacobi identity in ring-polymer phase space.

In this appendix, we provide the derivation of the main results of the paper, showing the connection and equivalence between the ring-polymer phase space and Hilbert space representation of time-correlation functions. A detailed step-by-step derivation can be found in Sec. IV of the supplementary material.

1. Auxiliary integrals

To demonstrate the connection between correlation functions in ring-polymer phase space and Hilbert space, it is convenient to introduce some auxiliary integrals that will simplify the algebra. To keep the presentation clear, we will use the shorthand notation

ÔW,jÔW(qj,pj)
(B1)

to denote the dependence on the jth coordinates of the (one-dimensional) Wigner–Weyl transform [Eq. (5)].29–32 Additionally, we will employ the operator identity

1̂=12πdqjdΔjdpjeipjΔjqj+Δj2qjΔj2.
(B2)

We introduce the “type-1” integral defined as

I1=ZN1(2π)NdqdpeβĤN̄1Nj=1NÔW,j,
(B3)

where Ô represents an arbitrary operator. Note that the integral defined above involves a ring-polymer phase-space average of a sum of one-dimensional Wigner transforms evaluated at particular phase-space points (qj, pj). However, recognizing the invariance of the Boltzmann operator [eβĤ]N̄ to cyclic permutations of the variables (i.e., invariance to the change qlql+1), the sum can be simplified to obtain an equivalent expression for I1 as follows:

I1=ZN1(2π)NdqdpeβĤN̄ÔW,N.
(B4)

It is straightforward to show that type-1 integrals of the form Eq. (B3) are connected to averages in Hilbert space. Specifically, introducing the definitions of the Wigner–Weyl transform and performing the integration over the ring-polymer phase-space coordinates (q, p) with the help of identity Eq. (B2), it follows that

I1=Z1TreβĤÔ.
(B5)

In other words, type-1 integrals correspond to averages in Hilbert space.

Similarly, we introduce the “type-2” integral defined as

I2=ZN1(2π)NdqdpeβĤN̄1Nj,k=1jkNÔW,jP̂W,k,
(B6)

where Ô and P̂ represent arbitrary operators. Note that the integral defined above involves a ring-polymer phase-space average of a sum of products of one-dimensional Wigner transforms evaluated at different phase-space points (qj, pj) and (qk, pk), with jk. However, recognizing the invariance of the Boltzmann operator [eβĤ]N̄ to cyclic permutations of the variables, the sum can be simplified to obtain an equivalent expression for I2, given by

I2=ZN1(2π)NdqdpeβĤN̄j=1N1ÔW,jP̂W,N.
(B7)

It is straightforward to perform the integration in Eq. (B7) over (q, p) with the help of identity Eq. (B2). In particular, the jth term in the sum is given by

Z1Tr(eβNĤ)jÔ(eβNĤ)NjP̂,
(B8)

from which it follows that

I2=j=1N1Z1Tr(eβNĤ)NjÔ(eβNĤ)jP̂.
(B9)

We remark that type-2 integrals correspond to (discrete) Kubo transform correlations between arbitrary operators Ô and P̂.45 

Similarly, we introduce the “type-3” integral defined as

I3=ZN1(2π)NdqdpeβĤN̄1Nj,k,l=1jklNÔW,jP̂W,kQ̂W,l,
(B10)

where Ô, P̂ and Q̂ represent arbitrary operators. Note that the integral defined above involves a ring-polymer phase-space average of a sum of products of one-dimensional Wigner transforms evaluated at different phase-space points (qj, pj), (qk, pk), and (ql, pl), with jk, jl, kl. However, recognizing the invariance of the Boltzmann operator [eβĤ]N̄ to cyclic permutations of the variables, the sum can be simplified to obtain an equivalent expression for I3,

I3=ZN1(2π)NdqdpeβĤN̄j,k=1jkN1ÔW,jP̂W,kQ̂W,N.
(B11)

Employing the identity Eq. (B2), integration over the ring-polymer phase-space coordinates can be performed. Two cases exist: For j < k, each term in the sum reduces to

Z1Tr(eβNĤ)jÔ(eβNĤ)kjP̂(eβNĤ)NkQ̂,
(B12)

whereas for j > k,

Z1Tr(eβNĤ)kP̂(eβNĤ)jkÔ(eβNĤ)NjQ̂.
(B13)

Combining the terms, it follows that

I3=j=1N1k=1j1Z1Tr(eβNĤ)NjÔ(eβNĤ)jkP̂(eβNĤ)kQ̂+j=1N1k=j+1N1Z1Tr(eβNĤ)NkP̂(eβNĤ)kjÔ(eβNĤ)jQ̂.
(B14)

We remark that type-3 integrals are related to (discrete) symmetrized Double Kubo transform correlations involving arbitrary operators Ô, P̂ and Q̂.45 

We remark that the procedure to define and evaluate “type-n” integrals In can be generalized to any order. However, as the order of the integral increases, the number of terms and possible permutations to be considered in the sums quickly grow and the derivation of the expressions becomes tedious. Nevertheless, clear trends can be identified: A type-n integral consists of ordered sums involving different imaginary time steps eβNĤ between all the possible permutations of operators [e.g., Eqs. (B9) and (B14)]. Such type-n integrals are discrete versions of symmetrized n-order Kubo transforms.

2. Sine correlation functions

To obtain the Hilbert space representation of the sine correlation functions [Eq. (21)], we used the definitions of generalized Wigner transforms [Eq. (4)] and sine coupling operator [Eq. (15)] to note that

Ô1NsÔ2N=i1Nj=1N[Ô1,Ô2]W,j.
(B15)

The sine coupling between two generalized Wigner transforms is just the sum of one-dimensional Wigner–Weyl transform of the commutator [Ô1,Ô2] evaluated at a particular ring-polymer coordinate (qj, pj).

Introducing the expression Eq. (B15) into Eq. (21), the sine correlation can be expressed as

Ô1sÔ2N=ZN1(2π)Ndqdp[eβĤ]N̄×i1Nj=1N[Ô1,Ô2]W,j.
(B16)

Besides an additional i/ factor, the sine correlation in Eq. (B16) represents a type-1 integral [Eq. (B3)] with Ô[Ô1,Ô2]. Therefore, it follows from Eq. (B5) that

Ô1sÔ2N=i[Ô1,Ô2],
(B17)

which proves Eq. (23). We remark that this equality holds for any finite N.

3. Cosine correlation functions

To obtain the Hilbert space representation of the cosine correlation functions [Eq. (22)], we used the definitions of generalized Wigner transforms [Eq. (4)] and cosine coupling operator [Eq. (16)] to note that

Ô1NcÔ2N=1N2j=1N12[Ô1,Ô2]+W,j+j,k=1jkN[Ô1]W,j[Ô2]W,k,,
(B18)

with [Ô1,Ô2]+=Ô1Ô2+Ô2Ô1 denoting an anti-commutator. The cosine coupling between two generalized Wigner transforms is composed of two terms: The first one is a sum of one-dimensional Wigner–Weyl transforms of the anti-commutator [Ô1,Ô2]+ evaluated at a particular bead coordinate; the second term, on the other hand, consists of sums of products between [Ô1]W and [Ô2]W evaluated at different bead coordinates.

Introducing the expression Eq. (B18) into Eq. (22), the cosine correlation can be decomposed into two terms as follows:

Ô1cÔ2N=12T1+T2,
(B19)

with

T1=ZN1(2π)Ndqdp[eβĤ]N̄1N2j=1N[Ô1,Ô2]+W,j,T2=ZN1(2π)Ndqdp[eβĤ]N̄1N2j,k=1jkN[Ô1]W,j[Ô2]W,k.
(B20)

Besides an additional 1/N factor, the T1 term in Eq. (B20) represents a type-1 integral [Eq. (B3)] with Ô[Ô1,Ô2]+, whereas the T2 term represents a type-2 integral [Eq. (B6)] with ÔÔ1 and P̂Ô2. Therefore, combining the terms, it follows from Eqs. (B5) and (B9) that

Ô1cÔ2N=Z1Nj=0NTr(eβNĤ)NjÔ1(eβNĤ)jÔ2,
(B21)

where the prime in ′ indicates that the first and last indices are weighted by one-half. Recognizing that in the N → ∞ limit the sum becomes a integral, it follows that

limNÔ1cÔ2N=1β0βdλÔ1(iλ)Ô2Ô1Ô2,
(B22)

which proves Eq. (24).

4. Sine–sine correlation functions

To obtain the Hilbert space representation of the sine–sine correlation functions [Eq. (27)], we used the definitions of generalized Wigner transforms [Eq. (4)] and sine coupling operator [Eq. (15)] and we note that

Ô1NsÔ2NsÔ3N=i21Nj=1N[[Ô1,Ô2],Ô3]W,j.
(B23)

The sine–sine coupling involves a sum of one-dimensional Wigner–Weyl transforms of the double commutator operator [[Ô1,Ô2],Ô3] evaluated at a particular bead coordinate.

Introducing the expression Eq. (B23) into Eq. (27), the sine–sine correlation can be expressed as

Ô1sÔ2sÔ3N=ZN1(2π)Ndqdp[eβĤ]N̄×i21Nj=1N[[Ô1,Ô2],Ô3]W,j.
(B24)

Besides an additional i/2 factor, the sine–sine correlation in Eq. (B24) represents a type-1 integral [Eq. (B3)] with Ô[[Ô1,Ô2],Ô3]. Therefore, it follows from Eq. (B5) that

Ô1sÔ2sÔ3N=i2[[Ô1,Ô2],Ô3],
(B25)

which proves Eq. (28). We remark that this equality holds for any finite N.

5. Sine–cosine correlation functions

To obtain the Hilbert space representation of the sine–cosine correlation functions [Eq. (29)], we note that

Ô1NsÔ2NcÔ3N=i1N2j=1N12[[Ô1,Ô2],Ô3]+W,j+j,k=1jkN[Ô1,Ô2]W,jÔ3W,k.
(B26)

The sine–cosine coupling is composed of two terms: the first one is a sum of one-dimensional Wigner–Weyl transforms evaluated at a particular bead coordinate of the symmetrized operator [Ô1,Ô2]Ô3; the second term, on the other hand, consists of sums of products between the commutator [[Ô1,Ô2]]W and [Ô3]W evaluated at different bead coordinates.

Introducing the expression Eq. (B26) into Eq. (29), the sine–cosine correlation can be decomposed into two terms as follows:

Ô1sÔ2cÔ3N=12T1+T2,
(B27)

with

T1=ZN1(2π)Ndqdp[eβĤ]N̄×i1N2j=1N[[Ô1,Ô2],Ô3]+W,j,
T2=ZN1(2π)Ndqdp[eβĤ]N̄×i1N2j,k=1jkN[Ô1,Ô2]W,j[Ô3]W,k.
(B28)

Besides an additional i/(Nℏ) factor, the T1 term in Eq. (B28) represents a type-1 integral [Eq. (B3)] with Ô[[Ô1,Ô2],Ô3]+, whereas the T2 term represents a type-2 integral [Eq. (B6)] with Ô[Ô1,Ô2] and P̂Ô3. Therefore, it follows from Eqs. (B5) and (B9) that

Ô1sÔ2cÔ3N=iZ1Nj=0NTr×(eβNĤ)Nj[Ô1,Ô2](eβNĤ)jÔ3,
(B29)

where the prime in ′ indicates that the first and last indices are weighted by one-half. Recognizing that in the N → ∞ limit the sum becomes a integral, it follows that

limNÔ1sÔ2cÔ3N=i1β0βdλe+λĤ[Ô1,Ô2]eλĤÔ3=i[Ô1,Ô2]Ô3,
(B30)

which proves Eq. (33).

6. Cosine–cosine correlation functions

To obtain the Hilbert space representation of the cosine–cosine correlation functions [Eq. (30)], we note that

Ô1NcÔ2NsÔ3N=1N3j=1N14[[Ô1,Ô2]+,Ô3]+W,j+j,l=1jlN12[Ô1,Ô2]+W,j[Ô3]W,l+j,k=1jkN12[Ô1]W,j[Ô2,Ô3]+W,k+j,k=1jkN12[Ô2]W,k[Ô1,Ô3]+W,j+j,k,l=1jk,jl,klN[Ô1]W,j[Ô2]W,k[Ô3]W,l.
(B31)

Note that the cosine–cosine coupling presents a more complex structure in which, depending on the sum term, different operators are evaluated at particular (equal or different) bead coordinates.

Introducing the expression Eq. (B31) into Eq. (30), the cosine–cosine correlation can be expressed as

Ô1cÔ2cÔ3N=14T1+12(T2+T2+T2)+T3,
(B32)

where each Tj term corresponds to a particular sum in Eq. (B31).

Besides a 1/N2 factor, all Tj terms in Eq. (B32) correspond to a type-1, type-2, or type-3 integral. For example, T1 corresponds to a type-1 integral with Ô[[Ô1,Ô2]+,Ô3]+, whereas the T2 terms correspond to type-2 integrals. The T3 term corresponds to a type-3 integral [Eq. (B10)] with ÔÔ1, P̂Ô2 and Q̂Ô3. Combining all the Tj terms together, it follows from Eqs. (B5), (B9) and (B14) that

Ô1cÔ2cÔ3N=Z1N2j=0Nk=0jTr(eβNĤ)NjÔ1×(eβNĤ)jkÔ2(eβNĤ)kÔ3+Z1N2j=0Nk=jN×Tr(eβNĤ)NkÔ2(eβNĤ)kjÔ1(eβNĤ)jÔ3,
(B33)

where the prime in ′ indicates that the first and last indices are weighted by one-half. Recognizing that the sums in the N → ∞ limit become an iterative double integral, it follows that

Ô1cÔ2cÔ3N=1β20βdλ00λ0dλ1Ô1(iλ0)Ô2(iλ1)Ô3+1β20βdλ0λ0βdλ1Ô2(iλ1)Ô1(iλ0)Ô3=1β20βdλ00βdλ1T̂βÔ1(iλ0)Ô2(iλ1)Ô3=Ô1Ô2Ô3,
(B34)

which proves Eq. (31). Note that we have introduced the imaginary time-ordering operator T̂β [Eq. (C5)] at the second equality to freely commute operators inside the integral sign.21,27

We define the Kubo transformed (KT) correlation function38 of two arbitrary operators as

Ô1Ô2=1β0βdλÔ1(iλ)Ô2,
(C1)

with Ô(iλ)=e+λĤÔeλĤ. It is straightforward to show that the KT presents the following symmetries:

Ô1Ô2=Ô2Ô1,
(C2)
Ô1Ô2*=Ô1Ô2.
(C3)

Note that the first symmetry allows expressing the KT in equivalent notations (i.e., the order of the argument is irrelevant), whereas the second symmetry indicates that the KT is a purely real correlation function if both Ô1 and Ô2 are Hermitian operators.

We define the (symmetrized) Double Kubo transformed (DKT) time correlation function as22,23,28,46

Ô1Ô2Ô3=1β20βdλ00βdλ1T̂βÔ1(iλ0)Ô2(iλ1)Ô3.
(C4)

Here, T̂β is an imaginary time-ordering operator that orders the product of operators so their imaginary time arguments increase from right to left, such that

T̂βÔ1(iλ0)Ô2(iλ1)=Ô1(iλ0)Ô2(iλ1)if λ0>λ1,Ô2(iλ1)Ô1(iλ0)if λ0<λ1.
(C5)

Note that T̂β ensures that there is no backward imaginary time propagation inside the integral.21 It is straightforward to show that the DKT satisfies the following symmetries:27 

Ô1Ô2Ô3=Ô2Ô1Ô3,
(C6)
Ô1Ô2Ô3=Ô2Ô3Ô1,
(C7)
Ô1Ô2Ô3*=Ô1Ô2Ô3.
(C8)

Note that the first and second symmetries indicate that the DKT is invariant to permutations and/or reordering of the arguments, allowing for different but equivalent notations. Moreover, the third symmetry implies that the DKT is a purely real correlation function for Hermitian operators Ô1, Ô2, and Ô3.

Similar considerations allow to extend the definition and notation to higher-order (fully symmetrized) Kubo transformed correlation functions.21 

Under the perturbative framework, the response of a system to a weak external perturbation can be described to nth order in terms of response functions defined as4,5

R(n)(t)=inTrÔn(tn)Ôn1(tn1),Ô1(t1),Ô0(t0),ρ̂.
(D1)

It is worth remarking that depending on the nature of the operators Ôj and the order n of the perturbation, the response functions can be related to transport coefficients, reaction rates, and several forms of linear and nonlinear vibrational spectroscopy (i.e., infrared, sum-frequency generation, two-dimensional Raman, two-dimensional terahertz-Raman).1,2,5

1. First-order

To first order (n = 1), the response function is given by

R(1)(t)=iB̂(t1),Â(t0).
(D2)

Using the Kubo identity38 

iÂ,ρ̂=ρ̂0βdλ0Â̇(iλ0),
(D3)

where Â̇(t)ddtÂ(t)=i[Ĥ,Â(t)] denotes a (real-)time derivative, Eq. (D2) can be also recast as

R(1)(t)=βddt0B̂(t1)Â(t0),
(D4)

in terms of the time derivative of a Kubo transformed TCF.

Equations (D2) and (D4) represent two alternative but equivalent ways of expressing the quantum first-order response functions in terms of a correlation function involving a commutator of operators or the time derivative of a Kubo time correlation function. We remark that these two correlations are mapped into the ring-polymer phase-space representation as the “sine” and “cosine” TCF, respectively.

2. Second-order

To second order (n = 2), the response function is given by

R(2)(t)=i2Ĉ(t2),B̂(t1),Â(t0),
(D5)

a time-correlation function involving a double commutator.

There are various alternative ways of recasting the second-order response function in terms of Kubo-like correlation functions.38 For example, using the Kubo identity [Eq. (D3)] Eq. (D5) can be recast as

R(2)(t)=βiddt0Ĉ(t2),B̂(t1)Â(t0),
(D6)

in terms of the time derivative of a Kubo-transform correlation function involving a commutator.

Alternatively, using the Double Kubo identity (see  Appendix E)

i2B̂,Â,ρ̂=ρ̂0βdλ00βdλ1T̂βÂ̇(iλ0)B̂̇(iλ1)iρ̂0βdλ0Â̇,B̂(iλ0),
(D7)

the response can be equivalently expressed as

R(2)(t)=β2d2dt0dt1Ĉ(t2)B̂(t1)Â(t0)+βiddt0Ĉ(t2)B̂(t1),Â(t0),
(D8)

where the second term is a time derivative of a KT correlation function [although different from the one presented in Eq. (D6)] and the first term represents the time derivatives of a (symmetrized) Double Kubo transformed time correlation function.

Equations (D5), (D6), and (D8) represent three equivalent ways of expressing the quantum second-order response functions. We remark that these correlations are mapped into the ring-polymer phase-space representation as the “sine–sine,” “sine-cosine,” and “cosine–cosine” TCF, respectively.

3. Third-order

To third-order (n = 3), the response function is given by

R(3)(t)=i3D̂(t3),Ĉ(t2),B̂(t1),Â(t0),
(D9)

a time-correlation function involving a triple commutator.

There are three alternative ways of recasting the third-order response function in terms of Kubo-like correlation functions. For example, using the Kubo identity [Eq. (D3)], Eq. (D9) can be recast as

R(3)(t)=βi2ddt0D̂(t3),Ĉ(t2),B̂(t1)Â(t0),
(D10)

in terms of the time derivative of a Kubo-transform correlation function involving a double commutator.

Alternatively, using the Double Kubo identity [Eq. (D7)] the response can be equivalently expressed as

R(3)(t)=β2id2dt0dt1D̂(t3),Ĉ(t2)B̂(t1)Â(t0)+βi2ddt0D̂(t3),Ĉ(t2)B̂(t1),Â(t0),
(D11)

where the first term involves time derivatives of a DKT involving a commutator and the second term represents a KT correlation function involving two commutators.

Finally, using the Triple Kubo identity (see  Appendix E)

i3Ĉ,B̂,Â,ρ̂=ρ̂0βdλ00βdλ10βdλ2T̂βÂ̇(iλ0)B̂̇(iλ1)Ĉ̇(iλ2)iρ̂0βdλ00βdλ1T̂βÂ̇(iλ0)B̂̇,Ĉ(iλ1)+B̂̇(iλ0)Â̇,Ĉ(iλ1)+Ĉ̇(iλ0)Â̇,B̂(iλ1)+i2ρ̂0βdλ0Â̇,B̂,Ĉ(iλ0),
(D12)

the response can be equivalently expressed as

R(3)(t)=β3d3dt0dt1dt2D̂(t3)Ĉ(t2)B̂(t1)Â(t0)+β2id2dt0dt1D̂(t3)Ĉ(t2),B̂(t1)Â(t0)+β2id2dt0dt1D̂(t3)Ĉ(t2),Â(t0)B̂(t1)+β2id2dt0dt2D̂(t3)B̂(t1),Â(t0)Ĉ(t2)+βi2ddt0D̂(t3)Ĉ(t2),B̂(t1),Â(t0).
(D13)

Equations (D9)(D11) and (D13) represent four equivalent ways of expressing the quantum third-order response functions. We remark that these correlations are mapped into the ring-polymer phase-space representation as different combinations of “sine” and “cosine” TCF [for example, a “sine–sine–sine” TCF for Eq. (D9) or a “sine–sine–cosine” TCF for Eq. (D10)].

In this appendix, we generalized the Kubo identity38 to higher orders. For simplicity in the derivation, we will use the notation

ÔλÔ(iλ)=eλĤÔeλĤ,
(E1)

to denote the dependence on imaginary time of an arbitrary operator Ô.

1. (Single) Kubo identity

We introduce the operator F̂A(λ), defined as

F̂A(λ)0λdλ0Â̇λ0,
(E2)

where a superdot on an operator denotes a time derivative, so

Â̇ddtÂ=i[Ĥ,Â].
(E3)

The integral over imaginary time in Eq. (E2) can be performed by noting that

Â̇λ=idÂλdλ,
(E4)

resulting in

F̂A(λ)=iÂλÂ.
(E5)

From Eq. (E5), it thus follows that

eβĤF̂A(β)=iÂ,eβĤ.
(E6)

Employing the definition of F̂A [Eq. (E2)], one obtains the (Single) Kubo Identity,38 

iÂ,eβĤ=eβĤ0βdλ0Â̇λ0.
(E7)

2. Double Kubo identity

We introduce the operator F̂AB(λ), defined as

F̂AB(λ)0λdλ00λdλ1T̂βÂ̇λ0B̂̇λ1,
(E8)

where T̂β is an imaginary time-ordering operator that orders the product of operators so their imaginary time arguments increase from right to left [Eq. (C5)].

Taking the derivative with respect to λ in Eq. (E8), it follows that

dF̂AB(λ)dλ=Â̇λF̂B(λ)+B̂̇λF̂A(λ),
(E9)

where we have used the time-ordering property of T̂β to order the product of operators and the definition of F̂O(λ) [Eq. (E2)].

Using the explicit expression Eq. (E5) for the operators F̂O(λ), it is straightforward to show that

F̂AB(β)=0βdλdF̂AB(λ)dλ=i0βdλÂ̇λB̂λ+B̂̇λÂλiF̂A(β)B̂+F̂B(β)Â.
(E10)

Recognizing that

iddλB̂λÂλ=B̂̇λÂλ+B̂λÂ̇λ,
(E11)

and, therefore, that

0βdλB̂̇λÂλ+B̂λÂ̇λ=iB̂βÂβB̂Â=B̂βF̂A(β)+F̂B(β)Â,
(E12)

where we have used Eq. (E5) to arrive at the last equality, Eq. (E10) can be written as

F̂AB(β)=i0βdλÂ̇λ,B̂λiF̂A(β)B̂B̂βF̂A(β).
(E13)

Recognizing that

eβĤF̂A(β)B̂B̂βF̂A(β)=iB̂,Â,eβĤ,
(E14)

and using the definition of F̂AB [Eq. (E8)], it follows from Eq. (E13) after reordering that

i2B̂,Â,eβĤ=eβĤ0βdλ00βdλ1T̂βÂ̇λ0B̂̇λ1ieβĤ0βdλÂ̇λ,B̂λ.
(E15)

This represents the Double Kubo Identity, a generalization of the Single Kubo identity to double commutators.

3. Triple Kubo identity

We introduce the operator F̂ABC(λ), defined as

F̂ABC(λ)0λdλ00λdλ10λdλ2T̂βÂ̇λ0B̂̇λ1Ĉ̇λ2,
(E16)

where T̂β the imaginary time-ordering operator.

Taking the derivative with respect to λ in Eq. (E16), it follows that

dF̂ABC(λ)dλ=Â̇λF̂BC(λ)+B̂̇λF̂AC(λ)+Ĉ̇λF̂AB(λ),
(E17)

where we have used the time-ordering property of T̂β to order the product of operators and the definition of F̂OO(λ) [Eq. (E8)].

Using the explicit expression Eq. (E13) for the operators F̂OO(λ), it is straightforward to show that

F̂ABC(β)=0βdλ0dF̂ABC(λ0)dλ0=i0βdλ00λ0dλ1Â̇λ0B̂̇λ1,Ĉλ1+B̂̇λ0Â̇λ1,Ĉλ1+Ĉ̇λ0Â̇λ1,B̂λ1i0βdλ0Â̇λ0F̂B(λ0)ĈĈλ0F̂B(λ0)+B̂̇λ0F̂A(λ0)ĈĈλ0F̂A(λ0)+Ĉ̇λ0F̂A(λ0)B̂B̂λ0F̂A(λ0).
(E18)

The expression in Eq. (E18) is expressed in terms of iterative double integrals of the form

0βdλ00λ0dλ1Ô̇λ0Ôλ1,
(E19)

where, for example, Ô=Â and Ô=B̂̇,Ĉ. Equivalent expressions can be obtained in terms of symmetrized double integrals of the form

0βdλ00βdλ1T̂βÔ̇λ0Ôλ1,
(E20)

by noting that

0βdλ00βdλ1T̂βÔ̇λ0Ôλ1=0βdλ00λ0dλ1Ô̇λ0Ôλ1+0βdλ0λ0βdλ1Ôλ1Ô̇λ0=0βdλ00λ0dλ1Ô̇λ0Ôλ1+0βdλ10λ1dλ0Ôλ1Ô̇λ0=0βdλ00λ0dλ1Ô̇λ0Ôλ1+0βdλ1Ôλ1F̂O(λ1).
(E21)

We have split the integral over λ1 and used the time-ordering property of the T̂β to obtain the first equality. At the second equality, we have interchanged the order of integration for the second term. The final equality follows from the definition of the operator F̂O in Eq. (E2).Applying Eq. (E21) to Eq. (E18) to express all iterative double integrals in terms of symmetrized double integrals, one obtains

F̂ABC(β)=i0βdλ00βdλ1T̂βÂ̇λ0B̂̇λ1,Ĉλ1+B̂̇λ0Â̇λ1,Ĉλ1+Ĉ̇λ0Â̇λ1,B̂λ1i0βdλ0Â̇λ0F̂B(λ0)ĈĈλ0F̂B(λ0)+B̂̇λ0F̂A(λ0)ĈĈλ0F̂A(λ0)+Ĉ̇λ0F̂A(λ0)B̂B̂λ0F̂A(λ0)+B̂̇λ0,Ĉλ0F̂A(λ0)+Â̇λ0,Ĉλ0F̂B(λ0)+Â̇λ0,B̂λ0F̂C(λ0)=i0βdλ00βdλ1T̂βÂ̇λ0B̂̇λ1,Ĉλ1+B̂̇λ0Â̇λ1,Ĉλ1+Ĉ̇λ0Â̇λ1,B̂λ1i20βdλ0Â̇λ0,B̂λ0,Ĉλ0Ĉ̇λ0B̂λ0Âλ0+Ĉλ0B̂̇λ0Âλ0+Ĉλ0B̂λ0Â̇λ0+B̂̇λ0Âλ0+B̂λ0Â̇λ0Ĉ+Ĉ̇λ0Âλ0+Ĉλ0Â̇λ0B̂+Ĉ̇λ0B̂λ0+Ĉλ0B̂̇λ0Â+i2F̂A(β)B̂Ĉ+F̂B(β)ÂĈ+F̂C(β)ÂB̂,
(E22)

where we have used the explicit expression for F̂O(λ0) [Eq. (E5)] to obtain the last equality.The previous expression can be simplified by using the operator identities

iddλ0Ôλ0Ôλ0=Ô̇λ0Ôλ0+Ôλ0Ô̇λ0,
(E23)

and

iddλ0Ôλ0Ôλ0Ôλ0=Ô̇λ0Ôλ0Ôλ0+Ôλ0Ô̇λ0Ôλ0+Ôλ0Ôλ0Ô̇λ0,
(E24)

resulting in

F̂ABC(β)=i0βdλ00βdλ1T̂βÂ̇λ0B̂̇λ1,Ĉλ1+B̂̇λ0Â̇λ1,Ĉλ1+Ĉ̇λ0Â̇λ1,B̂λ1i20βdλ0Â̇λ0,B̂λ0,Ĉλ0iddλ0Ĉλ0B̂λ0Âλ0+iddλ0B̂λ0Âλ0Ĉ+iddλ0Ĉλ0Âλ0B̂+iddλ0Ĉλ0B̂λ0Â+i2F̂A(β)B̂Ĉ+F̂B(β)ÂĈ+F̂C(β)ÂB̂=i0βdλ00βdλ1T̂βÂ̇λ0B̂̇λ1,Ĉλ1+B̂̇λ0Â̇λ1,Ĉλ1+Ĉ̇λ0Â̇λ1,B̂λ1i20βdλ0Â̇λ0,B̂λ0,Ĉλ0+i2F̂A(β)B̂ĈB̂βF̂A(β)Ĉ+ĈβB̂βF̂A(β)ĈβF̂A(β)B̂,
(E25)

where we have performed the integrals involving the ddλ0 terms at the last equality.

Recognizing that

eβĤF̂A(β)B̂ĈB̂βF̂A(β)Ĉ+ĈβB̂βF̂A(β)ĈβF̂A(β)B̂=iĈ,B̂,Â,eβĤ
(E26)

and using the definition of F̂ABC [Eq. (E16)], it follows from Eq. (E25) after reordering that

i3Ĉ,B̂,Â,eβĤ=eβĤ0βdλ00βdλ10βdλ2T̂βÂ̇λ0×B̂̇λ1Ĉ̇λ2ieβĤ0βdλ00βdλ1T̂βÂ̇λ0B̂̇λ1,Ĉλ1+B̂̇λ0Â̇λ1,Ĉλ1+Ĉ̇λ0Â̇λ1,B̂λ1+i2eβĤ0βdλ0Â̇λ0,B̂λ0,Ĉλ0.
(E27)

Equation (E27) represents the Triple Kubo Identity, a generalization of the Kubo identity to triple commutators.

1.
D.
Chandler
,
Introduction to Modern Statistical Mechanics
(
Oxford University Press
,
1987
).
2.
A.
Nitzan
,
Chemical Dynamics in Condensed Phases
(
Oxford University Press
,
2006
).
3.
D. A.
McQuarrie
,
Statistical Mechanics
(
University Science Books
,
2000
).
4.
S.
Mukamel
,
Principles of Nonlinear Optical Spectroscopy
(
Oxford University Press
,
1995
).
5.
M.
Cho
,
Two-Dimensional Optical Spectroscopy
(
CRC Press
,
2009
).
6.
J.
Maldacena
,
S. H.
Shenker
, and
D.
Stanford
,
J. High Energy Phys.
2016
,
106
.
7.
R. P.
Feynman
,
Rev. Mod. Phys.
20
,
367
(
1948
).
8.
R.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
McGraw-Hill Companies
,
1965
).
9.
D.
Chandler
and
P. G.
Wolynes
,
J. Chem. Phys.
74
,
4078
(
1981
).
10.
M.
Parrinello
and
A.
Rahman
,
J. Chem. Phys.
80
,
860
(
1984
).
11.
D. M.
Ceperley
,
Rev. Mod. Phys.
67
,
279
(
1995
).
12.
T. E.
Markland
and
M.
Ceriotti
,
Nat. Rev. Chem.
2
,
0109
(
2018
).
13.
N.
Makri
,
J. Math. Phys.
36
,
2430
(
1995
).
14.
M. S.
Litsarev
and
I. V.
Oseledets
,
J. Comput. Phys.
305
,
557
(
2016
).
15.
M. R.
Jèrgensen
and
F. A.
Pollock
,
Phys. Rev. Lett.
123
,
240602
(
2019
).
16.
A.
Bose
and
P. L.
Walters
,
J. Chem. Phys.
156
,
024101
(
2022
).
17.
T. J. H.
Hele
,
M. J.
Willatt
,
A.
Muolo
, and
S. C.
Althorpe
,
J. Chem. Phys.
142
,
134103
(
2015
).
18.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
100
,
5106
(
1994
).
19.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
101
,
6168
(
1994
).
20.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
Annu. Rev. Phys. Chem.
64
,
387
(
2013
).
21.
K. A.
Jung
,
P. E.
Videla
, and
V. S.
Batista
,
J. Chem. Phys.
151
,
034108
(
2019
).
22.
K. A.
Jung
,
P. E.
Videla
, and
V. S.
Batista
,
J. Chem. Phys.
153
,
124112
(
2020
).
23.
Z.
Tong
,
P. E.
Videla
,
K. A.
Jung
,
V. S.
Batista
, and
X.
Sun
,
J. Chem. Phys.
153
,
034117
(
2020
).
24.
T. J. H.
Hele
,
Mol. Phys.
115
,
1435
(
2017
).
25.
T. J. H.
Hele
,
M. J.
Willatt
,
A.
Muolo
, and
S. C.
Althorpe
,
J. Chem. Phys.
142
,
191101
(
2015
).
26.
T. J. H.
Hele
,
Mol. Phys.
114
,
1461
(
2016
).
27.
K. A.
Jung
,
P. E.
Videla
, and
V. S.
Batista
,
J. Chem. Phys.
148
,
244105
(
2018
).
28.
G.
Trenins
and
S. C.
Althorpe
,
J. Chem. Phys.
149
,
014102
(
2018
).
29.
M.
Hillery
,
R. F.
O’Connell
,
M. O.
Scully
, and
E. P.
Wigner
,
Phys. Rep.
106
,
121
(
1984
).
30.
H.-W.
Lee
,
Phys. Rep.
259
,
147
(
1995
).
31.
W. B.
Case
,
Am. J. Phys.
76
,
937
(
2008
).
32.
A.
Polkovnikov
,
Ann. Phys.
325
,
1790
(
2010
).
33.
S. C.
Althorpe
,
Eur. Phys. J. B
94
,
155
(
2021
).
34.
S.
Mukamel
,
Phys. Rev. E
68
,
021111
(
2003
).
35.
R.
Zwanzig
,
Nonequilibrium Statistical Mechanics
(
Oxford University Press
,
2001
).
36.
J. E.
Moyal
and
M. S.
Bartlett
,
Math. Proc. Cambridge Philos. Soc.
45
,
99
(
1949
).
37.
Q.
Shi
and
E.
Geva
,
J. Chem. Phys.
118
,
8173
(
2003
).
38.
R.
Kubo
,
J. Phys. Soc. Jpn.
12
,
570
(
1957
).
39.
D. R.
Reichman
,
P.-N.
Roy
,
S.
Jang
, and
G. A.
Voth
,
J. Chem. Phys.
113
,
919
(
2000
).
40.
S.
Mukamel
,
V.
Khidekel
, and
V.
Chernyak
,
Phys. Rev. E
53
,
R1
(
1996
).
41.
M.
Kryvohuz
and
J.
Cao
,
Phys. Rev. Lett.
96
,
030403
(
2006
).
42.
H.
Goldstein
,
C.
Poole
, and
J.
Safko
,
Classical Mechanics
(
Addison Weslet
,
2002
).
43.

In Ref. 21 (see also Refs. 17 and 28), we used a different notation to define the generalized Wigner–Weyl transform [Eq. (4)] applicable to position-only operators. In Sec. I D of the supplementary material, we show that both definitions are equivalent. The definition used here [Eq. (4)] is applicable to arbitrary operators of q̂ and p̂.

44.

For ring-polymer correlation functions involving only the sine coupling operator, the equality between Hilbert space and the ring-polymer phase-space representation holds for any value of N.

45.

Technically speaking, the end points of the Kubo integral discretization are not accounted for.

46.

We remark that in Ref. 27 (see also Refs. 21–23), we introduced different distinct DKT correlation functions and used a different notation for the symmetrized DKT given by Eq. (C4). In particular, the notation included a “sym” superscript. Because in the present paper, the only relevant DKT is the one defined by Eq. (C4), we drop the superscript in the notation.

Supplementary Material