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.

## I. INTRODUCTION

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 mechanics^{7,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-time^{17,18,24–26} and two-times^{22,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.

## II. IMAGINARY-TIME PATH-INTEGRAL PHASE SPACE

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 $H\u0302=p\u03022/2m+V(q\u0302)$. Nevertheless, the resulting expressions are also applicable to higher-dimensional systems.

### A. Ensemble averages

We focus on the thermal ensemble average of an operator $O\u0302=O(q\u0302,p\u0302)$,

where

is the partition function and *β* = 1/*k*_{B}*T*, with *T* the temperature and *k*_{B} the Boltzmann constant.

We introduce the ring-polymer phase-space average,^{17}

to compute expectation values in terms of phase-space averages with *∫d***q** ≡ *∫dq*_{1}⋯*∫dq*_{N}. Here, the *N* ring-polymer coordinates **q** = {*q*_{1}, …, *q*_{N}} are defined by the discretized imaginary time Feynman path subject to cyclic boundary conditions (i.e., *q*_{0} = *q*_{N}). Analogous definitions apply for the ring-polymer of momenta **p** = {*p*_{1}, …, *p*_{N}}.

The ring-polymer phase-space function $O\u0302N(q,p)$, introduced by Eq. (3), is the generalized Wigner–Weyl transform defined, as follows:^{43}

where

is the Wigner–Weyl transform of the operator $O\u0302$,^{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,

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

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\u2212\beta H\u0302N\u0304(q,p)$ defined, as follows:^{17,21,28}

with *β*_{N} = *β*/*N* and normalization constant given by

By evaluating the matrix elements $\u3008y|e\u2212\beta NH\u0302|x\u3009$ with the symmetric Trotter approximation in the *N* → ∞ limit, an explicit form for $e\u2212\beta H\u0302N\u0304$ can be obtained (see Sec. I C of the supplementary material) as follows:

with

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

with

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)

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 representation^{29–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}

### B. Quantum Liouvillian

Dynamical extensions of Eqs. (3) and (4) can be obtained for time-dependent operators, $O\u0302(t)=eiH\u0302t/\u210fO\u0302e\u2212iH\u0302t/\u210f$. In fact, the ring-polymer phase-space representation introduced in Sec. II A remains valid when introducing the substitution $O\u0302\u2192O\u0302(t)$ in Eq. (4) giving the time-dependent phase-space distribution $[O\u0302(t)]N(q,p)$. We remark that the Hamiltonian involved in the time evolution operator $e\xb1iH\u0302t/\u210f$ could differ from the one governing the statistics $e\u2212\beta H\u0302$, allowing the formalism introduced in this work to be applied to systems out of equilibrium.

The phase-space time evolution of $[O\u0302(t)]N(q,p)$ is obtained by extending the standard Liouvillian formalism^{34,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}

where the arrows indicate the direction in which the derivative is applied. Note that $\Lambda \u20e1N$ 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\u20e1$ and $c\u20e1$, which are defined as follows:

and

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

which represents a generalization of the Moyal expansion of the quantum Liouvillian^{29,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), namely^{21,24}

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

For systems at equilibrium, for which $e\u2212LNte\u2212\beta H\u0302N\u0304=e\u2212\beta H\u0302N\u0304$, 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 $[O\u0302]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 system^{17,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,

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 $[O\u0302(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.

## III. MULTI-TIME CORRELATION FUNCTIONS

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\u2212\beta H\u0302]N\u0304$ with the time evolution described by the ring-polymer propagator $eLNt[O\u0302]N$ to derive and introduce an exact ring-polymer phase-space representation of multipoint time correlation functions.

### A. Two-point time correlation functions

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

and

where $[O\u0302(t)]N=eLNt[O\u0302]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 $[A\u0302]N$ and $[B\u0302]N$ are evolved in time by the propagator $eLNt$ and coupled by a sine or cosine couplings operator.

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

with $[O\u03021,O\u03020]=O\u03021O\u03020\u2212O\u03020O\u03021$ 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

where we define the KT correlation function of two *arbitrary* operators as^{38}

(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 $A\u0302(t0)$ and $B\u0302(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}

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}

### B. Three-point time correlation functions

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:

Figure 4 shows a schematic representation of this correlation corresponding to the evolution of three observables $A\u0302$, $B\u0302$, and $C\u0302$ 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)].

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

so, the “sine–sine” correlation function is an exact path-integral representation of a correlation function involving a double commutation relation between observables $A\u0302(t0)$, $B\u0302(t1)$, and $C\u0302(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,

and

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

where we define the symmetrized DKT time correlation function as^{21,22,27,39}

with $T\u0302\beta $ 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,

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.

### C. General time correlation functions

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:

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

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:

where each $Jj\u20e1$, *j* = 1, …, *n*, represents either a sine coupling operator $(s\u20e1)$ or a cosine coupling operator $(c\u20e1)$. 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\u20e1$/$c\u20e1$), 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

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

and

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.

## IV. DISCUSSION

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.

### A. Invariance to imaginary time translation

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\u2212\beta H\u0302]N\u0304$ and Wigner–Weyl functions $[O\u0302]N$, as well as the Liouvillian $LN$ and Janus operator $\Lambda \u20e1N$, are invariant to the cyclic permutation *q*_{l} → *q*_{l+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 functions^{17,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.

### B. Symmetry

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\u20e1)$ in ring-polymer phase space, whereas the presence of a Kubo integral contributes with a cosine coupling operator $(c\u20e1)$. 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.

### C. Interference of quantum trajectories

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

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(\Lambda \u20e1N)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

with *x*_{i} = (*q*_{i}, *p*_{i}) 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 $Mij\u2026k(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.

### D. A word on notation

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\u2212\beta H\u0302]N\u0304c\u20e1[A\u0302(t0)]N$ or $[e\u2212\beta H\u0302]N\u0304s\u20e1[A\u0302(t0)]N$ factors instead. Note that if $[A\u0302(t0)]N$ is only first order in position or momenta, only the first term in the cosine/sine expansion will survive and $[e\u2212\beta H\u0302]N\u0304c\u20e1[A\u0302(t0)]N$ = $[e\u2212\beta H\u0302]N\u0304[A\u0302(t0)]N$, whereas $[e\u2212\beta H\u0302]N\u0304s\u20e1[A\u0302(t0)]N$ ∝ $[e\u2212\beta H\u0302]N\u0304\Lambda \u20e1N[A\u0302(t0)]N$. We remark that previously derived expressions for the KT and DKT by others^{17,28} and us^{21} 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.

### E. Emergence of classical limit

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

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

whereas $\u22c5cl$ 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\Lambda \u20e1O2\u2261\u2212{O1,O2}PB$ represents the negative of the classical Poisson bracket. In the case of three-point correlations, one obtains

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.

## V. CONCLUDING REMARKS

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 $\Lambda \u20e1N$ 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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

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

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

### APPENDIX A: PROPERTIES OF SINE AND COSINE COUPLING OPERATOR

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

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

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

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

### APPENDIX B: DERIVATION OF PATH-INTEGRAL PHASE-SPACE REPRESENTATION OF CORRELATION FUNCTIONS

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

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

We introduce the “type-1” integral defined as

where $O\u0302$ 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 (*q*_{j}, *p*_{j}). However, recognizing the invariance of the Boltzmann operator $[e\u2212\beta H\u0302]N\u0304$ to cyclic permutations of the variables (i.e., invariance to the change *q*_{l} → *q*_{l+1}), the sum can be simplified to obtain an equivalent expression for *I*_{1} as follows:

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

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

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

where $O\u0302$ and $P\u0302$ 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 (*q*_{j}, *p*_{j}) and (*q*_{k}, *p*_{k}), with *j* ≠ *k*. However, recognizing the invariance of the Boltzmann operator $[e\u2212\beta H\u0302]N\u0304$ to cyclic permutations of the variables, the sum can be simplified to obtain an equivalent expression for *I*_{2}, given by

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

from which it follows that

We remark that type-2 integrals correspond to (discrete) Kubo transform correlations between arbitrary operators $O\u0302$ and $P\u0302$.^{45}

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

where $O\u0302$, $P\u0302$ and $Q\u0302$ 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 (*q*_{j}, *p*_{j}), (*q*_{k}, *p*_{k}), and (*q*_{l}, *p*_{l}), with *j* ≠ *k*, *j* ≠ *l*, *k* ≠ *l*. However, recognizing the invariance of the Boltzmann operator $[e\u2212\beta H\u0302]N\u0304$ to cyclic permutations of the variables, the sum can be simplified to obtain an equivalent expression for *I*_{3},

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

whereas for *j* > *k*,

Combining the terms, it follows that

We remark that type-3 integrals are related to (discrete) symmetrized Double Kubo transform correlations involving arbitrary operators $O\u0302$, $P\u0302$ and $Q\u0302$.^{45}

We remark that the procedure to define and evaluate “type-n” integrals *I*_{n} 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\u2212\beta NH\u0302$ 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

The sine coupling between two generalized Wigner transforms is just the sum of one-dimensional Wigner–Weyl transform of the commutator $[O\u03021,O\u03022]$ evaluated at a particular ring-polymer coordinate (*q*_{j}, *p*_{j}).

Besides an additional *i*/*ℏ* factor, the sine correlation in Eq. (B16) represents a type-1 integral [Eq. (B3)] with $O\u0302\u2192[O\u03021,O\u03022]$. Therefore, it follows from Eq. (B5) that

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

with $[O\u03021,O\u03022]+=O\u03021O\u03022+O\u03022O\u03021$ 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 $[O\u03021,O\u03022]+$ evaluated at a particular bead coordinate; the second term, on the other hand, consists of sums of products between $[O\u03021]W$ and $[O\u03022]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:

with

Besides an additional 1/*N* factor, the *T*_{1} term in Eq. (B20) represents a type-1 integral [Eq. (B3)] with $O\u0302\u2192[O\u03021,O\u03022]+$, whereas the *T*_{2} term represents a type-2 integral [Eq. (B6)] with $O\u0302\u2192O\u03021$ and $P\u0302\u2192O\u03022$. Therefore, combining the terms, it follows from Eqs. (B5) and (B9) that

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

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

The sine–sine coupling involves a sum of one-dimensional Wigner–Weyl transforms of the double commutator operator $[[O\u03021,O\u03022],O\u03023]$ evaluated at a particular bead coordinate.

Besides an additional $i/\u210f2$ factor, the sine–sine correlation in Eq. (B24) represents a type-1 integral [Eq. (B3)] with $O\u0302\u2192[[O\u03021,O\u03022],O\u03023]$. Therefore, it follows from Eq. (B5) that

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

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 $[O\u03021,O\u03022]O\u03023$; the second term, on the other hand, consists of sums of products between the commutator $[[O\u03021,O\u03022]]W$ and $[O\u03023]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:

with

Besides an additional *i*/(*Nℏ*) factor, the *T*_{1} term in Eq. (B28) represents a type-1 integral [Eq. (B3)] with $O\u0302\u2192[[O\u03021,O\u03022],O\u03023]+$, whereas the *T*_{2} term represents a type-2 integral [Eq. (B6)] with $O\u0302\u2192[O\u03021,O\u03022]$ and $P\u0302\u2192O\u03023$. Therefore, it follows from Eqs. (B5) and (B9) that

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

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

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

where each *T*_{j} term corresponds to a particular sum in Eq. (B31).

Besides a 1/*N*^{2} factor, all *T*_{j} terms in Eq. (B32) correspond to a type-1, type-2, or type-3 integral. For example, *T*_{1} corresponds to a type-1 integral with $O\u0302\u2192[[O\u03021,O\u03022]+,O\u03023]+$, whereas the *T*_{2} terms correspond to type-2 integrals. The *T*_{3} term corresponds to a type-3 integral [Eq. (B10)] with $O\u0302\u2192O\u03021$, $P\u0302\u2192O\u03022$ and $Q\u0302\u2192O\u03023$. Combining all the *T*_{j} terms together, it follows from Eqs. (B5), (B9) and (B14) that

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

### APPENDIX C: KUBO TRANSFORM CORRELATION FUNCTIONS

We define the Kubo transformed (KT) correlation function^{38} of two *arbitrary* operators as

with $O\u0302(\u2212i\u210f\lambda )=e+\lambda H\u0302O\u0302e\u2212\lambda H\u0302$. It is straightforward to show that the KT presents the following symmetries:

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 $O\u03021$ and $O\u03022$ are Hermitian operators.

We define the (symmetrized) Double Kubo transformed (DKT) time correlation function as^{22,23,28,46}

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

Note that $T\u0302\beta $ 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}

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 $O\u03021$, $O\u03022$, and $O\u03023$.

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

### APPENDIX D: QUANTUM RESPONSE THEORY

Under the perturbative framework, the response of a system to a weak external perturbation can be described to *n*th order in terms of response functions defined as^{4,5}

It is worth remarking that depending on the nature of the operators $O\u0302j$ 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

Using the Kubo identity^{38}

where $A\u0302\u0307(t)\u2261ddtA\u0302(t)=i\u210f[H\u0302,A\u0302(t)]$ denotes a (real-)time derivative, Eq. (D2) can be also recast as

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

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

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

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

the response can be equivalently expressed as

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.

#### 3. Third-order

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

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

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

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)

the response can be equivalently expressed as

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

### APPENDIX E: GENERALIZED KUBO IDENTITIES

In this appendix, we generalized the Kubo identity^{38} to higher orders. For simplicity in the derivation, we will use the notation

to denote the dependence on imaginary time of an arbitrary operator $O\u0302$.

#### 1. (Single) Kubo identity

We introduce the operator $F\u0302A(\lambda )$, defined as

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

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

resulting in

From Eq. (E5), it thus follows that

#### 2. Double Kubo identity

We introduce the operator $F\u0302AB(\lambda )$, defined as

where $T\u0302\beta $ 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

where we have used the time-ordering property of $T\u0302\beta $ to order the product of operators and the definition of $F\u0302O(\lambda )$ [Eq. (E2)].

Using the explicit expression Eq. (E5) for the operators $F\u0302O(\lambda )$, it is straightforward to show that

Recognizing that

and, therefore, that

Recognizing that

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\u0302ABC(\lambda )$, defined as

where $T\u0302\beta $ the imaginary time-ordering operator.

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

where we have used the time-ordering property of $T\u0302\beta $ to order the product of operators and the definition of $F\u0302OO\u2032(\lambda )$ [Eq. (E8)].

Using the explicit expression Eq. (E13) for the operators $F\u0302OO\u2032(\lambda )$, it is straightforward to show that

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

where, for example, $O\u0302=A\u0302$ and $O\u0302\u2032=B\u0302\u0307,C\u0302$. Equivalent expressions can be obtained in terms of *symmetrized* double integrals of the form

by noting that

We have split the integral over *λ*_{1} and used the time-ordering property of the $T\u0302\beta $ 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\u0302O$ in Eq. (E2).Applying Eq. (E21) to Eq. (E18) to express all iterative double integrals in terms of symmetrized double integrals, one obtains

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

and

resulting in

where we have performed the integrals involving the $dd\lambda 0$ terms at the last equality.

Recognizing that

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

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

## REFERENCES

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\u0302$ and $p\u0302$.

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

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

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.