The simulation of multidimensional vibrational spectroscopy of condensed-phase systems including nuclear quantum effects is challenging since full quantum-mechanical calculations are still intractable for large systems comprising many degrees of freedom. Here, we apply the recently developed double Kubo transform (DKT) methodology in combination with ring-polymer molecular dynamics (RPMD) for evaluating multi-time correlation functions [K. A. Jung et al., J. Chem. Phys. 148, 244105 (2018)], providing a practical method for incorporating nuclear quantum effects in nonlinear spectroscopy of condensed-phase systems. We showcase the DKT approach in the simulation of the fifth-order two-dimensional (2D) Raman spectroscopy of Lennard-Jones liquids as a prototypical example, which involves nontrivial nonlinear spectroscopic observables of systems described by anharmonic potentials. Our results show that the DKT can faithfully reproduce the 2D Raman response of liquid xenon at high temperatures, where the system behaves classically. In contrast, liquid neon at low temperatures exhibits moderate but discernible nuclear quantum effects in the 2D Raman response compared to the responses obtained with classical molecular dynamics approaches. Thus, the DKT formalism in combination with RPMD simulations enables simulations of multidimensional optical spectroscopy of condensed-phase systems that partially account for nuclear quantum effects.

Ultrafast spectroscopy has been demonstrated to be powerful in accessing dynamical information of chromophores in liquids, such as time-dependent fluorescence Stokes shift for probing solvation dynamics1–4 and two-dimensional (2D) vibrational spectra for correlating vibrations at different spatial and temporal points.5–7 Most ultrafast spectroscopic techniques share a common feature as in the traditional pump–probe approach: a time-dependent laser pulse pumps the system either vibrationally or electronically, and at a delayed time, the probe pulse is used to read out the molecular response to the disturbance caused by the pump pulse. At first glance, this pump–probe process involves an explicitly time-dependent Hamiltonian, but if the pulses are weak enough, one can express the nonequilibrium optical response as an equilibrium time-correlation function (TCF) according to the widely applicable linear response theory.1,2 In this context, the fluctuation-dissipation theorem8 of linear response theory connects the equilibrium fluctuations with the nonequilibrium relaxation of observables and can be applied to treat multiple pumps as in multidimensional spectroscopy.9 Consequently, the nonlinear spectroscopic signal is expressed in terms of the convolution between the optical response function and the laser-pulse envelopes in the time domain. The response function carries all the information of the physical system and provides dynamical information such as coupling between vibrational states and energy transfer pathways in light-harvesting systems.10,11

Like most of the condensed-phase measurements, linear spectroscopy—such as electronic, IR, THz, Raman, and sum-frequency generation spectra—constitutes a dimensionality reduction from 3N to 2 for an N-atom system, whose response function is usually an ensemble average of some spectroscopic observables that reflects highly averaged molecular signatures. When we go to multidimensional spectroscopy, e.g., two-dimensional spectra,6,12 extra structural and dynamical information can be revealed by the correlation between multiple dimensions. A famous example is the coupling between a pair of amide I local modes in proteins via hydrogen-bonding, resulting in cross peaks in the 2D IR spectra.13,14

The fifth-order 2D Raman spectroscopy15–23 introduces an interaction (Hint = −E1ΠE2/2, where E1/2 are ultrafast electric fields and Π is the many-body polarizability) prior to the conventional third-order optical Kerr effect (OKE) spectroscopy,24–27 where the nonlinear order corresponds to the number of perturbing electric fields, some of which may come from the same laser pulse. The 2D Raman pulse sequence is shown in Fig. 1. Although the 2D Raman signal is relatively weak due to the high nonlinear order compared to the third-order cascaded contamination,18 it has been reported for CS2,19,20,28,29 benzene,30 and formamide liquids,31 as well as photoactive proteins.23 Recently, the single-beam spectrally controlled 2D Raman technique has been developed to overcome the cascading problem, leading to a rebirth of interest in 2D Raman spectroscopy.32 Theoretical studies of 2D Raman spectroscopy offer unique capabilities for studying vibrational mode coupling,33,34 anharmonicity of interactions,35–37 and dephasing processes in liquids.38,39

FIG. 1.

Pulse sequence of the fifth-order two-dimensional (2D) Raman spectroscopy. The first pair of Raman pump electric fields E1, E2 is applied at time 0, and the second pair of Raman pump electric fields E3, E4 is applied at time t1; then, the fifth-order Raman signal Es is scattered off from the impeding probe electric field E5 at time t1 + t2. Including the signal field, the fifth-order Raman scattering is a six-wave-mixing nonlinear optical process.

FIG. 1.

Pulse sequence of the fifth-order two-dimensional (2D) Raman spectroscopy. The first pair of Raman pump electric fields E1, E2 is applied at time 0, and the second pair of Raman pump electric fields E3, E4 is applied at time t1; then, the fifth-order Raman signal Es is scattered off from the impeding probe electric field E5 at time t1 + t2. Including the signal field, the fifth-order Raman scattering is a six-wave-mixing nonlinear optical process.

Close modal

The simulation of 2D Raman spectroscopy has been challenging due to the multidimensional characteristics of the TCF involved. In terms of classical mechanics, the 2D Raman response function involves two nested Poisson brackets, although the outmost Poisson bracket can be eliminated if the phase-space density is the equilibrium Boltzmann distribution.35,40 The computation of the remaining Poisson bracket is challenging since it is a quantity related to the stability of classical trajectories in phase space and is usually hard to converge for intrinsically chaotic many-body molecular systems such as liquids.41–43 Thus, the direct evaluation of the Poisson bracket or equivalently the stability matrix is computationally difficult. Ma and Stratt systematically investigated the fifth-order 2D Raman spectra using all-atom classical molecular dynamics (MD) simulation and the instantaneous-normal-mode (INM) theory 44 and elucidated the sensitivity of the 2D Raman signal with respect to dynamical anharmonicity and nonlinear polarizability.17,35,36,45,46 Additionally, mode-coupling theory seems capable of reproducing some features in 2D Raman spectroscopy of simple liquids.47 Tanimura and co-workers recently developed a hybrid equilibrium–nonequilibrium molecular dynamics approach to address the high computational overhead, where nonequilibrium molecular dynamics simulations are used to evaluate the Poisson bracket by mimicking the physical effect of the perturbing electric fields.48–50 It is also worth noting that DeVane, Ridley, Space, and Keyes (DRSK) proposed an alternative way to cast the classical 2D response function into ordinary TCFs based on the harmonic approximation,51–54 as discussed later in this paper.

The classical representation in atomistic simulations of liquids may introduce errors in the spectrum due to the neglect of nuclear quantum effects (NQEs)55 such as tunneling and zero-point energy, which could be significant in systems containing light particles and at low temperatures. The NQEs of condensed matter are typically modeled using Feynman’s path-integral formulation of quantum mechanics. For time-independent or thermodynamic properties such as structures or free energies, quantum statistical mechanics is isomorphic to sampling in an extended configuration space, treating each physical particle as a P-bead (or P-replica) ring polymer where the adjacent beads are connected with harmonic springs and the beads with the same bead index (1, …, P) belonging to different physical particles, interacting with each other according to the external potential, forming P “parallel bead universes.” When P, we will have exact quantum statistics, but, in practice, we choose a finite number P with convergence. The sampling in this extended ring-polymer space could be implemented fully classically using Monte Carlo as in path-integral Monte Carlo (PIMC) or molecular dynamics with additional fictitious bead masses and momenta as in path-integral molecular dynamics (PIMD).56–58 Since eβĤ with inverse temperature β = 1/kBT (kB is Boltzmann’s constant) can be considered as a propagator eiĤt/ with imaginary time t = −iβℏ, PIMC and PIMD are commonly addressed as imaginary-time path-integral approaches, which are only used to sample the extended configuration space to average time-independent quantum properties. The averages converge quickly since the imaginary-time path integral does not have oscillatory amplitudes as in the real-time propagator.

Approximate methods such as centroid molecular dynamics (CMD)59–61 and ring-polymer molecular dynamics (RPMD)62–64 have been proposed to incorporate NQE in time-dependent properties such as TCFs, where the trajectories of the ring polymers are viewed as real-time dynamics with delocalized nuclei. With recently developed thermostat and propagation algorithms, CMD and RPMD have shown great potential for the simulation of linear spectroscopy in liquids.65–68 However, a general approach for incorporating NQE into multidimensional spectroscopy and multi-time correlation functions was introduced only recently, as shown for the 2D spectroscopic response function expressed in terms of double Kubo transformed (DKT) correlation functions, estimated by using RPMD.69 

In this paper, we focus on the simulation of 2D Raman spectroscopy of atomic liquids using RPMD. Using the DKT approach, the 2D response function in the frequency domain can be broken into real (symmetric) and imaginary (asymmetric) parts of the DKT correlation functions, where the symmetric part can be directly simulated via RPMD and is assumed to be dominant. Neglecting the asymmetric part is the only approximation in the DKT approach and was tested to be valid in 1-dimensional model systems.69 In this paper, we apply and test the DKT approach to more realistic and anharmonic Lennard-Jones (LJ) liquids, including Xe and Ne, to account for NQEs in their 2D Raman spectra obtained from RPMD simulations and to test the validity of the approximation involved in the DKT approach for atomic liquids. Lennard-Jones potentials are simple atomic liquid models; however, the anharmonic interaction in 3N-dimensional atomistic space and the nonlinear many-body polarizability observable make these calculations rather challenging.

This paper is organized as follows: Sec. II summarizes the theory of 2D Raman spectroscopy evaluation of DKT correlation functions with RPMD as well as classical MD including the equilibrium approach involving the stability matrix, the DRSK TCF approach, and the hybrid equilibrium–nonequilibrium approach. Section III describes the Lennard-Jones liquid models (Xe and Ne) and the simulation details. Section IV discusses the comparison between the classical and quantum 2D Raman spectra for Xe and Ne. Concluding remarks are given in Sec. V. More detailed accounts of some of the mathematical derivations are given in the  Appendix.

The optical observable in fifth-order 2D Raman spectroscopy can be expressed in terms of the two-time response function as9 

R(5)t1,t2=i2TrΠ^t1+t2Π^t1,Π^(0),ρ^eq.
(1)

Here, Π^ represents the many-body polarizability tensor, Â(t) = eiĤt/ÂeiĤt/ is the operator evolved in the Heisenberg picture, [Â,B^]=ÂB^B^Â is the quantum commutator, and ρ^eq=eβĤ/Z is the equilibrium density operator with system Hamiltonian Ĥ, inverse temperature β = 1/kBT, and partition function Z = Tr[eβĤ]. By selecting the polarization of the perturbing electric fields Ei coupled to the Raman interactions (Fig. 1), different elements of the many-body polarizability tensor Π can be probed. In this study, we focus on the fully polarized 2D Raman response for which the element Π^=Π^=Π^zz in Eq. (1) (generalization to other polarization conditions is straightforward).

A classical mechanical limit of the 2D response function can be obtained by invoking the mapping relation between the quantum commutator and the Poisson bracket,70 

iÂ,B^A,B,
(2)

where

A,B=ArBpBrAp.
(3)

Here, A = A(r, p) and B = B(r, p) represent classical dynamical variables, r = (r1, …, rN) are the 3N-dimensional coordinates, and p = (p1, …, pN) are the corresponding momenta of an N-atom system. Using this quantum–classical correspondence and noting that {A, ρeq} = −βȦρeq with Ȧ = dA/dt,71 the 2D Raman response function can then be expressed as35 

R(5)t1,t2=Πt1+t2,Πt1,Π(0)=βΠ̇(t1)Π0,Πt2,
(4)

where ⟨·⟩ = ∫ drdpρeq represents a classical ensemble average with Hamiltonian H, phase-space density ρeq = eβH/Z, and classical partition function Z = ∫ drdpeβH.

We use three different approaches to evaluate the classical response function in Eq. (4). The first approach is based on equilibrium MD simulations and relies on the fact that the observables depend only on positions r, i.e., Π^=Π^(r). The Poisson bracket in the classical 2D response function can be expressed as

Π0,Πt2=Π(r0)r0Π(rt2)rt2Jrp(t2),
(5)

and the classical 2D response is thus

R(5)t1,t2=βΠ̇(t1)Π(r0)r0Π(rt2)rt2Jrp(t2).
(6)

Here, Jrp(t) = rt/p0 is the stability matrix (or Jacobian matrix) satisfying the equations of motion J̈rp(t)=(1/m)D(t)Jrp(t) with initial condition Jrp(0) = 0, mass m, and dynamical matrix (or instantaneous Hessian) D(t)=V(r)rt.35 The stability matrix describes the sensitivity of positions at time t, rt, if the initial momenta, p0, are disturbed. In the equilibrium MD approach, the stability matrix has to be converged along with the spectroscopic observables in Eq. (6). In practice, the stability matrix is difficult to evaluate due to the fact that the classical trajectories of many-body systems are intrinsically chaotic.41 Fortunately, this calculation is feasible for simple liquids and is considered as the exact classical result for the proposal of this study.35 

The second approach that we consider is the hybrid equilibrium–nonequilibrium MD,48–50 where the Poisson bracket is estimated by averaging pairs of positively and negatively perturbed nonequilibrium MD trajectories, namely,

{Π(0),Π(t2)}=1E1E2ΔtΠ+Π(0)(t2)ΠΠ(0)(t2).
(7)

Here, Π±Π(0) (t2) represents the many-body polarizability at time t2 calculated from nonequilibrium trajectories that are generated by an ultrafast Raman perturbation Hint = ±(−E1ΠE2/2) at time t = 0, with external electric fields E1/2 of time duration Δt.40,72 The extra forces due to the Raman perturbation are given by ΔF = −rHint = ±E1rΠE2/2. The 2D Raman response function in the hybrid equilibrium–nonequilibrium is then given by

R(5)t1,t2=βE1E2ΔtΠ̇(t1)Π+Π(0)(t2)ΠΠ(0)(t2),
(8)

where Π̇(t1) is evaluated from an equilibrium trajectory. This approach is computationally more efficient than the equilibrium approach since it does not involve the convergence problem of the stability matrix and provides virtually the same results (vide infra).

The third approach is based on the approximate methodology developed by DeVane, Ridley, Space, and Keyes (DRSK) to compute 2D Raman spectroscopy in the classical limit.51 Starting with the exact quantum 2D response function given by Eq. (1) in Ref. 51, it was shown that invoking a harmonic reference potential and using the high temperature (classical) limit, the 2D response function can be approximated as

RDRSK(5)(t1,t2)=β22C(t1,t2)t1t2122C(t1,t2)t12,
(9)

where C(t1, t2) represents the classical two-time correlation function of the fluctuations of the many-body polarizability, namely,

C(t1,t2)=δΠ(0)δΠt1δΠt1+t2,
(10)

with δΠ=ΠΠ. We remark that the response given by Eq. (1) or (4) is invariant under the replacement Π → δΠ; however, Eq. (9) only holds for polarizability fluctuations.51 The response in Eq. (9) represents an approximation of the Poisson bracket contribution by the classical response function54 and with the advantage that no explicit perturbations need to be simulated and only equilibrium trajectories (that do not involve the stability matrix) are necessary to converge the TCF, providing a tractable theory of the classical nonlinear response.51–54 

A quantum-mechanically exact expression for the 2D response function [Eq. (1)] can be derived in terms of the Double Kubo Transform (DKT) time correlation function. In Ref. 69, it was shown that the 2D response can be recast in the Fourier domain as

R̃(5)(ω1,ω2)=Q+(ω1,ω2)K̃sym(ω1,ω2)+Q(ω1,ω2)K̃asym(ω1,ω2),
(11)

where the Q factors are defined as

Q±(ω1,ω2)=12Q1(ω1,ω2)±Q1(ω1,ω2),
(12)

with

Q1(ω1,ω2)=β2eβω11ω1ω2(ω1ω2)eβω1ω2eβω2ω1+ω1ω2
(13)

and K̃sym(ω1,ω2) and K̃asym(ω1,ω2) are the double Fourier transforms of the symmetric (real) and asymmetric (imaginary) DKT correlation functions, namely,

Ksym(t1,t2)=2ReΠ;Π(t1);Π(t1+t2),
(14)
Kasym(t1,t2)=2iImΠ;Π(t1);Π(t1+t2).
(15)

Here, the DKT correlation function is defined as69,73

Π;Π(t1);Π(t1+t2)=1Zβ20βdλ0λdλTreβĤΠ^(iλ)×Π^(t1iλ)Π^(t1+t2).
(16)

It is noted that in Ref. 69, we used a different choice of time variables t1, t2 in the response function. However, with a change of variables (t1, t2) → (t1, t1 + t2) and (ω1, ω2) → (ω1ω2, ω2), the expressions derived in Ref. 69 can be adapted to the current time definition. Detailed derivations using the relative times are provided in the supplementary material.

Equation (11) is quantum-mechanically exact. However, approximations need to be made for practical applications to realistic condensed-phase systems. Reference 69 demonstrated that the main contribution to the response is given by the symmetric part of the DKT, at least for systems with harmonic or mildly anharmonic potentials. We therefore neglect the asymmetric term in Eq. (11) and approximate the 2D Raman response, in what we call the DKT approach, as69 

R̃DKT(5)(ω1,ω2)Q+(ω1,ω2)K̃sym(ω1,ω2).
(17)

Moreover, semiclassical methodologies including multi-time Matsubara dynamics74,75 and multi-time RPMD62,69 were developed to evaluate the symmetric DKT correlation function Ksym(t1, t2).69,75 In this study, we use RPMD to approximate the symmetric DKT in terms of polarizability fluctuations as

Ksym(t1,t2)=limPδΠ(0)δΠ(t1)δΠ(t1+t2)P,
(18)

where

δΠ(0)δΠ(t)δΠ(t)P=1ZP(2π)fdR0dP0eβPHP(R0,P0)×δΠP(R0)δΠP(Rt)δΠP(Rt).
(19)

Here, f = 3NP is the total number of degrees of freedom with N as the number of particles and P as the number of replicas used to discretize the Boltzmann distribution, βP = β/P, and R=({rj(α)})(j=1,,N;α=1,,P) are the ring-polymer coordinates (with a similar definition for the conjugate momenta P). ZP is the P-bead path-integral representation of the partition function,

ZP=1(2π)fdRdPeβPHP(R,P).
(20)

AP is the ring-polymer representation of any observable A,

AP(R)=1Pα=1PA(r1(α),,rN(α)).
(21)

HP is the ring-polymer Hamiltonian,

HP(R,P)=HP0(R,P)+α=1PV(r(α)),
(22)

where V(r) is the external potential energy, while the free ring-polymer Hamiltonian is

HP0(R,P)=j=1Nα=1P(pj(α))22mj+12mjωP2rj(α)rj(α1)2,
(23)

with ωP = 1/(βP) and rj(0)rj(P). In Eq. (19), Rt denotes the ring-polymer coordinates at time t evolved from R0 at time t = 0 using the Hamiltonian HP.

The RPMD approximation [Eq. (18)] can be shown to be an exact path-integral discretization of the symmetric DKT in the limit t1, t2 → 0 and to be exact for all times for linear operators in a harmonic potential. Moreover, it can be derived from the more general multi-time Matsubara dynamics.76 The advantage of the DKT formulation with respect to the classical formulations presented in Sec. II A is that the TCF in Eq. (19) is sampled from the correct quantum Boltzmann distribution, hence providing a tractable theory that includes NQE in the calculation of 2D response functions.69 

The atomic liquid systems are described by the pairwise Lennard-Jones (LJ) potential, 77 

V(r)=i<jNu(rij)=i<jN4εσrij12σrij6,
(24)

where rij = |rij| = |rirj| is the interatomic distance between atom i and atom j, ε is the energy depth parameter, and σ is the length parameter of the LJ interaction. We study two model systems, including (i) a LJ fluid representing xenon with parameters ε/kB = 222 K, σ = 4.099 Å, and mass m = 131.293 amu (with reduced time unit τLJ=mσ2/ε=3.47 ps) at the phase point T*TkB/ε = 1.00, ρ*ρσ3 = 0.80 and (ii) a LJ fluid representing neon with parameters ε/kB = 35.6 K, σ = 2.749 Å, and m = 20.1797 amu (τLJ = 2.27 ps) at the phase point T* = 0.84, ρ* = 0.78. Both phase points correspond to the liquid phase and are chosen to represent a liquid at high and low temperatures, respectively. It should be pointed out that for classical LJ systems at a given phase point, all equilibrium and dynamical properties are independent of the atomic species, provided that reduced units are used. However, when nuclear quantum effects are included, the results depend on Planck’s constant (in reduced units, */mϵσ2) that is dependent on the mass of the particle (* ≈ 0.009 95 for Xe and * ≈ 0.094 53 for Ne).78 

Our MD simulations of liquid Xe included N = 108 atoms at temperature T* = 1.00 and density ρ* = 0.80. RPMD simulations of liquid Xe including N = 108 atoms under the same conditions represented each atom by P = 16 beads. The simulations are performed within a cubic periodic box of length 5.130σ using the velocity Verlet integrator79 with a time step of δt = 0.0005τLJ = 1.73 fs. Equilibration is achieved in a 1.73 ns NVE simulation, as checked with standard protocols.80 Liquid configurations are typically sampled every five time steps. In the equilibrium classical MD simulation, 1 × 108 liquid configurations from 500 independent trajectories of 1.73 ns were averaged. In the hybrid equilibrium–nonequilibrium MD simulation, 1 × 108 liquid configurations generated from equilibrium MD were sampled to generate 1 × 108 pairs of positively and negatively perturbed nonequilibrium trajectories by applying extra forces with external electric fields E1 = E2 = 5.0 V/Å. In the DRSK TCF calculation, 1 × 108 liquid configurations from 500 independent trajectories of 1.73 ns were averaged to obtain δΠ(0)δΠt1δΠt1+t2. RPMD simulations were implemented in the framework of thermostated RPMD (TRPMD),65 where a Langevin thermostat with parameter γ(k) = 2ωk was applied to the normal modes with positive ring-polymer frequencies {ωk > 0} (i.e., no thermostat applied to the centroid with zero ring-polymer frequency to preserve the Hamiltonian dynamics), and 4 × 108 ring-polymer configurations from 200 independent trajectories of 17.3 ns were averaged to compute the TCF. In all cases, the TCF was computed for |t1,2| < 1.73 ps (while we only plotted the time range of −1.0 ps to 1.0 ps in Sec. IV).

Simulations of liquid Ne were performed analogously with N = 108, T* = 0.84, ρ* = 0.78 for classical MD simulations and N = 64 represented by P = 64 beads for RPMD, in a cubic periodic box of length 5.173σ, integrated using the velocity Verlet algorithm with a time step of δt = 0.001τLJ = 2.27 fs. Equilibration was achieved in a 2.27 ns NVE simulation. A total of 1 × 108 ring-polymer configurations from 500 independent trajectories (2.27 ns each) were averaged to obtain δΠ(0)δΠt1δΠt1+t2P. The TCF was computed with a maximum t1 and t2 of 2.04 ps. Table I summarizes the simulation parameters.

TABLE I.

The simulation parameters of liquid Xe (at T* = 1.00, ρ* = 0.80) and liquid Ne (at T* = 0.84, ρ* = 0.78) including the system size (N), number of beads per atom (P), time step (δt), number of trajectories and configurations averaged, and kinetic energy per particle per kBT.

NPδt/τLJTraj.Config.βKE⟩/N
Xe (RPMD) 108 16 0.0005 200 4 × 108 1.50 ± 0.03 
Xe (MD) 108 0.0005 500 1 × 108 1.50 ± 0.07 
Ne (RPMD) 64 64 0.001 500 1 × 108 1.85 ± 0.04 
Ne (MD) 108 0.001 500 1 × 108 1.50 ± 0.08 
NPδt/τLJTraj.Config.βKE⟩/N
Xe (RPMD) 108 16 0.0005 200 4 × 108 1.50 ± 0.03 
Xe (MD) 108 0.0005 500 1 × 108 1.50 ± 0.07 
Ne (RPMD) 64 64 0.001 500 1 × 108 1.85 ± 0.04 
Ne (MD) 108 0.001 500 1 × 108 1.50 ± 0.08 

The many-body polarizability Π needed for the calculation of the fifth-order full-polarized 2D Raman response R(5) (t1, t2) was computed using the dipole-induced-dipole (DID) model,81 where the total system polarizability Π is expressed as the sum of effective site polarizabilities πj (j = 1, …, N),

Π=j=1Nπj.
(25)

The effective site polarizability is expressed in terms of N gas-phase (or isolated) atomic polarizabilities αj (j = 1, …, N) and interaction-induced terms involving the dipole–dipole interaction tensor T,

πj=αj1+kjNT(rjk)πk
(26)

with

T(r)=1r=3r^r^1r3,
(27)

where r^=r/r,r=|r|,and (rjk=rkrj).

In the case of LJ atomic liquids, each atomic polarizability tensor αj = αj1 is a constant diagonal 3 × 3 matrix. The atomic polarizability of Xe αXe = 4.11 Å3 and Ne αNe = 0.396 Å3 was taken from Ref. 82. Note that the interaction-induced terms involve the dipole–dipole tensor T(r′), which is a nonlinear function of the interatomic distances, thus providing a sensitive probe of liquid structure rearrangements.83 

The effective site polarizability [Eq. (26)] can be computed through self-consistent iterations starting with the initial guess πj = αj or by the more computational costly matrix inversion. To reduce the computational overhead, we limit ourselves to the first-order DID approximation

πj=αj+αjkjNT(rjk)αk.
(28)

Previous work showed that the first-order DID approximation can capture the polarizability dynamics of this atomic liquid model, although it underestimates the absolute magnitudes of the response functions.84,85

To further reduce the computational costs of converging multi-time correlation functions, we took advantage of the fact that in isotropic systems, such as liquids, averages of products of tensors can be expressed in terms of rotational invariants, namely,86 

AzzBzzCzz¯=1105Tr(A)Tr(B)Tr(C)+8105TP(A,B,C)+2105Tr(A)PP(B,C)+Tr(B)PP(A,C)+Tr(C)PP(A,B),
(29)

where the trace, pair product, and triple product are, respectively, defined as

Tr(A)=iAii,
(30)
PP(A,B)=i,jAijBij,
(31)
TP(A,B,C)=i,j,kAijBikCjk.
(32)

In our case, Azz, Bzz, and Czz correspond to polarizability fluctuations evaluated at different times [see Eqs. (10) and (18)]. It should be mentioned that in the hybrid equilibrium–nonequilibrium approach, the isotropic symmetry of the liquid is broken due to the explicit perturbation from the external electric fields applied along the z direction; thus, rotational averages cannot be applied for this method.40 

It is worth noting that varying polarization conditions that correspond to different tensor elements in the 2D Raman response could be used to identify different coupling mechanisms by selecting out specific molecular degrees of freedom,45,46,87 and the corresponding rotational average expressions do not depend on whether the dynamics is treated classically or quantum-mechanically.86 However, polarization selectivity is out of the scope of this work.

We begin by assessing the significance of nuclear quantum effects in LJ liquids. To this end, in Table I, we present the average kinetic energy per particle for liquid Xe and liquid Ne at the corresponding state points. The average kinetic energy is computed using the PIMD virial estimator given by88 

KEP=3N2β+12Pj=1Nα=1Prj(α)r¯jVr1(α),,rN(α)rj(α)P,
(33)

where r¯j=1Pα=1Prj(α) is the centroid position of atom j. The first term in Eq. (33) represents the kinetic energy associated with the motion of the ring-polymer centroids r¯j, whereas the second term accounts for delocalization and quantum fluctuations of the ring-polymer beads with respect to their centroids and could be substantial in systems with significant NQEs. Note that in the classical P = 1 case, the kinetic energy estimator reduces to the classical kinetic energy, namely, ⟨KE1 = 3N/(2β). Thus, values of βKEP/N that deviate from 1.5 can be used to assess the importance of NQE. As can be seen from Table I, for the case of liquid Xe at ρ* = 0.80 and T* = 1.0, there is no difference between the quantum kinetic energy and the classical result, indicating that Xe behaves as a classical liquid. On the other hand, the quantum average kinetic energy for liquid Ne at ρ* = 0.84 and T* = 0.78 presents an increase of ∼25% with respect to the classical result, indicating significant NQE for liquid Ne at low temperatures.

The influence of NQE on the structural properties of liquids can be inferred from the radial distribution function (RDF). In Figs. 2(a) and 2(b), we present g(r) computed from PIMD and classical MD simulations for both Ne and Xe liquids. As can be observed, the RDF for liquid xenon is almost identical, showing no influence of NQE in the structure of the liquid. On the other hand, noticeable differences can be seen in the peak height and position of the first solvation shell in liquid Ne, indicating a moderate influence of NQE on the coordination environment around Ne atoms.

FIG. 2.

Radial distribution functions obtained from ring-polymer molecular dynamics (RPMD, red solid curve) and molecular dynamics (MD, black dashed curve) for (a) liquid Xe (at T* = 1.00, ρ* = 0.80) and (b) liquid Ne (at T* = 0.84, ρ* = 0.78). The velocity–velocity time correlation functions obtained from RPMD (red solid curve) and MD (black dashed curve) for (c) liquid Xe and (d) liquid Ne. The number of atoms, N, and number of beads per atom, P, are indicated. Lennard-Jones reduced units (σ and τLJ) are used. Results are obtained by averaging over 107 configurations for RPMD of Ne, and 108 configurations for the other cases.

FIG. 2.

Radial distribution functions obtained from ring-polymer molecular dynamics (RPMD, red solid curve) and molecular dynamics (MD, black dashed curve) for (a) liquid Xe (at T* = 1.00, ρ* = 0.80) and (b) liquid Ne (at T* = 0.84, ρ* = 0.78). The velocity–velocity time correlation functions obtained from RPMD (red solid curve) and MD (black dashed curve) for (c) liquid Xe and (d) liquid Ne. The number of atoms, N, and number of beads per atom, P, are indicated. Lennard-Jones reduced units (σ and τLJ) are used. Results are obtained by averaging over 107 configurations for RPMD of Ne, and 108 configurations for the other cases.

Close modal

As a final descriptor of NQE on LJ liquids, we performed RPMD simulations to compute the velocity autocorrelation function defined as64,89

Cvv(t)=1Nj=1Nv¯j(0)v¯j(t),
(34)

which represents an approximation to the Kubo transformed velocity autocorrelation function. In the previous equation, v¯j=1Pα=1Pvj(α) corresponds to the centroid velocity of atom j. By comparing the RPMD velocity autocorrelation function with its classical P = 1 counterpart, it is possible to assess the influence of NQE in the dynamics of LJ liquids. In Figs. 2(c) and 2(d), we present the results for Cvv(t) computed with RPMD and MD simulations for both Ne and Xe. For the case of Xe, both correlation functions are almost identical, showing no noticeable NQE. On the contrary, for liquid Ne, the quantum TCF decays ∼0.05τLJ = 0.1 ps slower than the classical one, indicating that NQEs have a moderate influence on the dynamics of the liquid.

The overall conclusion of the previous analysis is that liquid Xe at ρ* = 0.80 and T* = 1.0 does not present any noticeable NQE and can be described as a classical liquid. This fluid thus provides a benchmark system for which classical dynamics should be accurate and serves as a test of the performance of the quantum dynamical method in the classical limit. On the other hand, liquid Ne at ρ* = 0.84 and T* = 0.78 exhibits moderate NQE on the liquid structure and the dynamics and hence serves as a test model for assessing NQE in nonlinear spectroscopies. It is important to remark that despite some influence of NQE being seen in liquid Ne, they are still not significant enough to induce exchange effects of nuclei like in superfluid 4He droplets below 1 K.57 

Having established the influence of NQE on LJ liquids, we now analyze the performance of the DKT approach to simulate 2D Raman spectroscopy. As an initial test, we considered the fully polarized 2D Raman response of (classic-like) liquid Xe at T* = 1.0 and ρ* = 0.80. Since the DKT approach is based on the TCF of fluctuations of the polarizability, in Fig. 3, we present contour plots of δΠ(0)δΠt1δΠt1+t2P computed with RPMD and classical (P = 1) MD simulations. The TCF smoothly decays to zero in a timescale of ∼1 ps from its initial value, showing little structure after ∼0.6 ps. Note that the central symmetry of the TCF [i.e., f(t1, t2) = f(−t1, −t2)] 69 is evident, which can be traced back to the time reversal symmetry. The comparison of the TCF obtained with RPMD and MD shows no discernible difference, which is consistent with the results obtained from the velocity autocorrelation function (Fig. 2).

FIG. 3.

Time correlation functions of many-body polarizability fluctuation δΠ(0)δΠt1δΠt1+t2 for liquid Xe obtained with (a) RPMD and (b) MD as well as for liquid Ne obtained with (c) RPMD and (d) MD. TCF are normalized.

FIG. 3.

Time correlation functions of many-body polarizability fluctuation δΠ(0)δΠt1δΠt1+t2 for liquid Xe obtained with (a) RPMD and (b) MD as well as for liquid Ne obtained with (c) RPMD and (d) MD. TCF are normalized.

Close modal

In the left panels of Fig. 4, we present the time-domain 2D Raman response functions for liquid Xe computed using the DKT approach [panel (a)]. We also include results of the classical response computed with the DRSK approximation [panel (b)], the stability matrix approach [panel (c)], and the hybrid equilibrium–nonequilibrium approach [panel (d)]. A comparison between the bottom two panels in Fig. 4 shows that the stability matrix and the hybrid equilibrium–nonequilibrium 2D spectra are identical, providing a confirmation that Poisson brackets evaluated by nonequilibrium MD are also numerically exact on the classical level.40 The classical time-domain 2D Raman response is predominantly characterized by a negative signal center at t1 ∼ 50 fs, t2 ∼ 340 fs, in consistency with previous results.35,36 Note also the lack of an echo signal along the diagonal t1 = t235,36,45,46 and that R(5) (t1, t2) = 0 along the t2 = 0 axis.

FIG. 4.

Fifth-order 2D Raman response functions for liquid Xe (left panels) and liquid Ne (right panels) computed via the DKT approach with RPMD [(a) and (e)], classical DRSK TCF approximation with MD [(b) and (f)], equilibrium classical MD approach involving stability matrix [(c) and (g)], and hybrid equilibrium–nonequilibrium MD [(d) and (h)]. Response functions are normalized.

FIG. 4.

Fifth-order 2D Raman response functions for liquid Xe (left panels) and liquid Ne (right panels) computed via the DKT approach with RPMD [(a) and (e)], classical DRSK TCF approximation with MD [(b) and (f)], equilibrium classical MD approach involving stability matrix [(c) and (g)], and hybrid equilibrium–nonequilibrium MD [(d) and (h)]. Response functions are normalized.

Close modal

The comparison of the DKT and the DRSK response with the classical result shows that the characteristic features of the 2D Raman response are fairly well reproduced by these approximated approaches.52 In particular, the overall lineshape and decay timescales of the spectra are in good agreement, with a negative signal centered at t1 ∼ 80 fs, t2 ∼ 340 fs, and a response that vanishes (within numerical error) along the t2 = 0 axis. However, both DKT and DRSK responses show a faster decay along the t2 axis for t2 < 0.4 ps, and a different spectral decay in the region t1 < 0.5 ps, t2 < 0.2 ps, somehow overestimating the intensity of the signal in this region (for better visualization of this effect, in Fig. 5, we present a 3-dimensional plot with 2D contour projection of the DKT and classical MD 2D Raman spectra of liquid Xe). This latter result is akin to the one found in one-dimensional systems and can be attributed to the neglect of the asymmetric part in the DKT approximation.69 

FIG. 5.

Fifth-order 2D Raman response functions for liquid Xe calculated via the (a) DKT approach with RPMD and (b) equilibrium MD approach. Response functions are normalized.

FIG. 5.

Fifth-order 2D Raman response functions for liquid Xe calculated via the (a) DKT approach with RPMD and (b) equilibrium MD approach. Response functions are normalized.

Close modal

To better quantify the differences in liquid Xe 2D Raman spectra obtained by different approaches, in Fig. 6, we present time slices of the time-domain 2D Raman responses calculated via DKT with RPMD, DRSK with MD, and equilibrium MD. The time slice at t2 = 0.34 ps with error bars of the DKT response and DRSK response [Fig. 6(g)] is shown in Fig. 9(a). Both approximate approaches show a small time shift in the position of the peak of the negative feature with respect to the converged classical response,52 an effect that is more evident for time slices at short-time [Figs. 6(a) and 6(b)]. The overestimation of the initial positive peak is also evident for times t2 < 0.2 ps. At longer times, the decay of the 2D response is accurately reproduced by both DKT and DRSK approaches. Given the overall performance, the results presented here confirm that both DKT and DRSK approaches provide computationally tractable theories for the calculation of 2D response functions in condensed-phase systems.

FIG. 6.

Time slices along t2 with t1 fixed [(a)–(d)], along t1 with t2 fixed [(e)–(g)], and the diagonal one with t1 = t2 (h) of the fifth-order 2D Raman response functions for liquid Xe calculated via the equilibrium MD approach (black), DRSK TCF approximation (cyan), and DKT with RPMD (red). The sampling rate is Δt = 5δt = 8.65 fs.

FIG. 6.

Time slices along t2 with t1 fixed [(a)–(d)], along t1 with t2 fixed [(e)–(g)], and the diagonal one with t1 = t2 (h) of the fifth-order 2D Raman response functions for liquid Xe calculated via the equilibrium MD approach (black), DRSK TCF approximation (cyan), and DKT with RPMD (red). The sampling rate is Δt = 5δt = 8.65 fs.

Close modal

It is interesting to note that the response obtained using the DKT approach and the classical DRSK approximation for liquid Xe is quite similar, even though both methods do not appear to be equivalent. To rationalize the connection between these two methodologies, it is instructive to note that the state point considered for liquid Xe corresponds to a high temperature (β → 0). In the high temperature limit, the intra-polymer harmonic spring forces in the RPMD Hamiltonian [Eq. (23)] become so large that the ring polymer collapses to a single point and the multi-time RPMD TCF effectively reduces to a multi-time classical correlation function, namely, Ksym(t1, t2) → C(t1, t2).69 Moreover, it is straightforward to show that in this high temperature limit (see the  Appendix),

limβ0Q+(ω1,ω2)=β23ω1(ω12ω2).
(35)

Therefore, in the high temperature limit, the symmetrized DKT approximation to the response is given by

limβ0R̃DKT(ω1,ω2)=2β2312ω12ω1ω2C̃(ω1,ω2),
(36)

which can be analytically Fourier transformed into the time domain to obtain

limβ0RDKT(5)(t1,t2)=2β232C(t1,t2)t1t2122C(t1,t2)t12.
(37)

Note that other than a constant multiplicative factor (which is immaterial in the normalized response shown in Fig. 6), Eq. (37) is equivalent to the DRSK approximation introduced by Eq. (9), demonstrating that both approaches are equivalent (to a scaling factor) in the high temperature limit. Note that the novel relation established in Eq. (37) is general for any two-time TCF and does not depend on the system under consideration, providing one of the main findings of this paper. Moreover, since the DKT approach stems from a rigorous quantum-mechanical expression valid at any temperature, unlike the DRSK approach, which is derived in the classical (high temperature) limit with harmonic potential assumptions, the DKT approach provides a clear advantage over the DRSK method, allowing one to include NQE in the computation of 2D response functions.

As a final test of the performance of the DKT approach, we analyzed the 2D Raman response of liquid Ne at T* = 0.84 and ρ* = 0.78, which presents moderate NQE. RPMD and classical results for the two-time correlation function of the fluctuations of the polarizability δΠ(0)δΠt1δΠt1+t2 are presented in the right panels of Fig. 3. The TCFs present similar characteristics to the ones obtained for liquid Xe, albeit a faster decay, which can be correlated with the shorter characteristic timescale τLJ for Ne as compared with Xe. The comparison between RPMD and MD TCFs shows a small difference in the lineshape of the correlation, with a slower decay for the former correlations, akin to the NQE observed in the velocity autocorrelation function (Fig. 2).

Figure 4 (right panels) shows contour plots of the time-domain 2D Raman response functions for liquid Ne obtained using different methodologies. Similar to liquid Xe, the response is characterized by a main negative peak at around t1 = 50 fs, t2 = 230 fs. However, the decay of the signals is faster in Ne than in Xe, in accordance with τLJ(Ne) < τLJ(Xe). In terms of the approximate approaches, both DKT and DRSK methods capture the overall lineshape and decay timescale of the response, although the decay along t2 is faster and the intensity of the short-time positive signal is again overestimated (see Fig. 7 for a 3-dimensional plot of the DKT and classical MD 2D Raman spectra of liquid Xe).

FIG. 7.

Fifth-order 2D Raman response functions for liquid Ne calculated via the (a) DKT approach with RPMD and (b) equilibrium MD approach. Response functions are normalized.

FIG. 7.

Fifth-order 2D Raman response functions for liquid Ne calculated via the (a) DKT approach with RPMD and (b) equilibrium MD approach. Response functions are normalized.

Close modal

Time slices of the 2D response of liquid Ne are presented in Fig. 8. Similar features as in the case of liquid Xe can be observed for the approximate methodologies in comparison with the classical result, namely, a small shift in the position of the negative peak and an overestimation of the intensity of the positive peak at shorter times. However, note that in this case, the response obtained with DKT and DRSK methods is not equivalent as evident by a slightly faster decay of the quantum response [Figs. 8(d) and 8(f)–8(h)]. Figure 9(b) shows the time slice at t2 = 0.23 ps with error bars of the DKT response and DRSK response [Fig. 8(g)], revealing that the difference between them is statistically distinguishable. Given the connection between both approaches in the high temperature limit (see Sec. IV B), the subtle difference between the DKT and DRSK response can be correlated with the different information encoded in the two-time RPMD TCF in contrast to the classical TCF (see Fig. 3). Since RPMD can encode NQE in multi-time correlation functions, as shown in Ref. 69, for simple model systems for which exact quantum calculations can be done, we therefore ascribed the faster decay in the DKT response function to NQE. Given the difficulty of performing exact quantum calculations for condensed-phase systems, our results show the DKT approach as an effective methodology to incorporate NQE in 2D spectroscopy.

FIG. 8.

Time slices along t2 with t1 fixed [(a)–(d)], along t1 with t2 fixed [(e)–(g)], and the diagonal one with t1 = t2 (h) of the fifth-order 2D Raman response functions for liquid Ne [Figs. 4(e)–4(g)] calculated via equilibrium the MD approach (black), DRSK TCF approximation (cyan), and DKT with RPMD (red). Here, the sampling rate is Δt = 5δt = 11.35 fs.

FIG. 8.

Time slices along t2 with t1 fixed [(a)–(d)], along t1 with t2 fixed [(e)–(g)], and the diagonal one with t1 = t2 (h) of the fifth-order 2D Raman response functions for liquid Ne [Figs. 4(e)–4(g)] calculated via equilibrium the MD approach (black), DRSK TCF approximation (cyan), and DKT with RPMD (red). Here, the sampling rate is Δt = 5δt = 11.35 fs.

Close modal
FIG. 9.

Time slices with error bars of the fifth-order 2D Raman response functions calculated via the DRSK TCF approximation (cyan) and DKT with RPMD (red) in liquid Xe (a) and Ne (b), corresponding to Figs. 6(g) and 8(g).

FIG. 9.

Time slices with error bars of the fifth-order 2D Raman response functions calculated via the DRSK TCF approximation (cyan) and DKT with RPMD (red) in liquid Xe (a) and Ne (b), corresponding to Figs. 6(g) and 8(g).

Close modal

In the present work, we have shown that the recently developed double Kubo transform (DKT) approach69 in combination with multi-time ring-polymer molecular dynamics (RPMD) provides a practical numerical method for the incorporation of nuclear quantum effects in the simulation of fifth-order 2D Raman spectroscopy of atomic liquids. Two different systems were studied. For liquid Xe at T* = 1.00, ρ* = 0.80, where NQEs are not important and the system behaves classically, the DKT 2D Raman response is comparable with the response obtained with classical methodologies. Moreover, in this high temperature regime, the DKT approach reduces (besides a constant multiplicative factor) to the classical method developed by DeVane et al.51 for the simulation of 2D Raman spectroscopy, providing a general result that is independent of the system and is one of the main findings of this paper. On the other hand, the results for liquid Ne at low temperatures (T* = 0.84, ρ* = 0.78) calculated using the DKT approach reveal moderate but discernible nuclear quantum effects in the 2D Raman response when compared with classical simulations. Given the fact that the 2D Raman spectroscopy depends on the many-body polarizability, which is a nonlinear operator of the nuclear coordinates, our results demonstrate the feasibility of the DKT approach to incorporate nuclear quantum effects in the simulation of nonlinear spectroscopy for condensed-phase systems.

Our study was focused on the 2D Raman response function of LJ atomic liquids. Although NQEs are observed in the 2D response for liquid Ne at low temperatures, the effects are rather small. It would be interesting to compare the calculated 2D Raman spectra with experimental measurements, which is currently unavailable. It would also be interesting to analyze how NQEs influence the spectroscopic response in systems that might show more pronounced quantum effects, such as para-H2,89 glassy systems,90,91 water,92 molecular liquids,93 or organic photovoltaic systems.94 Moreover, the DKT approach can be readily applied to other types of vibrational spectroscopy whose response can be expressed as a multi-time correlation involving three operators, such as 2D Raman-THz,49,95,96 where NQE can be more pronounced.97 Future work will explore the performance of the DKT formalism for these systems.

Finally, we would like to remark possible improvements to the current DKT formalism, which involves two approximations: (1) the neglect of the asymmetric part in the exact DKT expansion of the 2D response function and (2) the use of RPMD to estimate the symmetrized DKT correlation function. The development of practical methods to incorporate the asymmetric DKT will extend the usage of the theory for the simulation of systems with significant anharmonicity. On the other hand, the dynamics generated by RPMD could be replaced by more accurate quantum dynamical methods. We have recently shown that the symmetrized DKT could be approximated by a multi-time version of Matsubara dynamics.75 Although not a practical method for condensed-phase systems, Matsubara dynamics provides a benchmark starting point for the future development of Boltzmann preserving semiclassical approximations to multi-time correlation functions. Additionally, extensions of other accurate quantum dynamical methods such as linearized path-integral methods98–100 could be used to give rise to better estimation of multi-time correlation functions. Finally, extension of the DKT framework for the simulation of electronic spectroscopies, such as 2D solute-pump/solvent-probe,40,84,85,93 or other vibrational spectroscopies, such as 2D THz101 or 2D IR spectroscopy,6,9,13,14 are in principle possible. Work in these directions is underway and will be reported in future publications.

See the supplementary material for the derivation of the DKT formalism using the relative times.

The data that support the findings of this study are available within the article and its supplementary material.

X.S. and P.E.V. thank Francesco Paesani for helpful discussions. X.S. acknowledges support from NYU Shanghai, the National Natural Science Foundation of China (Grant No. 21903054), the Shanghai Sailing Program (Grant No. 19YF1435600), and the Program for Eastern Young Scholar at Shanghai Institutions of Higher Learning. Computing resources were provided by the NYU Shanghai High Performance Computing Center. V.S.B. acknowledges support from the NSF under Grant No. CHE-1900160 and generous allocation of computer time from the Yale High Performance Computing Center.

The high temperature limit of the Q± factors presented in Eq. (12) can be obtained by performing a Taylor expansion in β of the exponentials in Q1 [Eq. (13)],

Q1(ω1,ω2)=22βω1(βω1)2/2+O[(β)3]1β(ω1+ω2)/3+O[(β)2]=22βω112(βω1)2+O[(β)3]×1+13β(ω1+ω2)+O[(β)2]=22βω116(β)2ω1(ω12ω2)+O[(β)3].
(A1)

Noting the definition of Q± [Eq. (12)], the high temperature limit of the Q factors is given by

limβ0Q+(ω1,ω2)β23ω1(ω12ω2),
(A2)
limβ0Q(ω1,ω2)2βω1.
(A3)
1.
R. M.
Stratt
and
M.
Maroncelli
,
J. Phys. Chem.
100
,
12981
(
1996
).
2.
E. A.
Carter
and
J. T.
Hynes
,
J. Chem. Phys.
94
,
5961
(
1991
).
3.
R.
Jimenez
,
G. R.
Fleming
,
P. V.
Kumar
, and
M.
Maroncelli
,
Nature
369
,
471
(
1994
).
4.
G. R.
Fleming
and
M.
Cho
,
Annu. Rev. Phys. Chem.
47
,
109
(
1996
).
5.
P.
Hamm
and
M.
Zanni
,
Concepts and Methods of 2D Infrared Spectroscopy
(
Cambridge University Press
,
Cambridge
,
2011
).
6.
M.
Cho
,
Chem. Rev.
108
,
1331
(
2008
).
7.
A. T.
Krummel
,
P.
Mukherjee
, and
M. T.
Zanni
,
J. Phys. Chem. B
107
,
9165
(
2003
).
8.
R.
Kubo
,
Rep. Prog. Phys.
29
,
255
(
2002
).
9.
S.
Mukamel
,
Principles of Nonlinear Optical Spectroscopy
(
Oxford University Press
,
New York
,
1995
).
10.
T.
Brixner
,
J.
Stenger
,
H. M.
Vaswani
,
M.
Cho
,
R. E.
Blankenship
, and
G. R.
Fleming
,
Nature
434
,
625
(
2005
).
11.
J.
Li
,
H.
Qian
,
H.
Chen
,
Z.
Zhao
,
K.
Yuan
,
G.
Chen
,
A.
Miranda
,
X.
Guo
,
Y.
Chen
,
N.
Zheng
,
M. S.
Wong
, and
J.
Zheng
,
Nat. Commun.
7
,
10749
(
2016
).
12.
T. l. C.
Jansen
,
S.
Saito
,
J.
Jeon
, and
M.
Cho
,
J. Chem. Phys.
150
,
100901
(
2019
).
13.
M.
Reppert
and
A.
Tokmakoff
,
Annu. Rev. Phys. Chem.
67
,
359
(
2016
).
14.
Z.
Ganim
,
H. S.
Chung
,
A. W.
Smith
,
L. P.
DeFlores
,
K. C.
Jones
, and
A.
Tokmakoff
,
Acc. Chem. Res.
41
,
432
441
(
2008
).
15.
R. F.
Loring
and
S.
Mukamel
,
J. Chem. Phys.
83
,
2116
(
1985
).
16.
Y.
Tanimura
and
S.
Mukamel
,
J. Chem. Phys.
99
,
9496
(
1993
).
17.
A.
Ma
and
R. M.
Stratt
,
Phys. Rev. Lett.
85
,
1004
(
2000
).
18.
D. A.
Blank
,
L. J.
Kaufman
, and
G. R.
Fleming
,
J. Chem. Phys.
111
,
3105
(
1999
).
19.
D. A.
Blank
,
L. J.
Kaufman
, and
G. R.
Fleming
,
J. Chem. Phys.
113
,
771
(
2000
).
20.
O.
Golonzka
,
N.
Demirdöven
,
M.
Khalil
, and
A.
Tokmakoff
,
J. Chem. Phys.
113
,
9893
(
2000
).
21.
D. A.
Blank
,
G. R.
Fleming
,
M.
Cho
, and
A.
Tokmakoff
, in
Ultrafast Infrared and Raman Spectroscopy
, edited by
M. D.
Fayer
(
Marcel Dekker
,
New York
,
2001
).
22.
K. J.
Kubarych
,
C. J.
Milne
, and
R. J. D.
Miller
,
Int. Rev. Phys. Chem.
22
,
497
(
2003
).
23.
H.
Kuramochi
,
S.
Takeuchi
,
H.
Kamikubo
,
M.
Kataoka
, and
T.
Tahara
,
Sci. Adv.
5
,
eaau4490
(
2019
).
24.
N. T.
Hunt
,
A. A.
Jaye
, and
S. R.
Meech
,
Phys. Chem. Chem. Phys.
9
,
2167
(
2007
).
25.
Q.
Zhong
and
J. T.
Fourkas
,
J. Phys. Chem. B
112
,
15529
(
2008
).
26.
S.
Ryu
and
R. M.
Stratt
,
J. Phys. Chem. B
108
,
6782
(
2004
).
27.
H. E.
Bailey
,
Y.-L.
Wang
, and
M. D.
Fayer
,
J. Chem. Phys.
149
,
044501
(
2018
).
28.
L. J.
Kaufman
,
J.
Heo
,
L. D.
Ziegler
, and
G. R.
Fleming
,
Phys. Rev. Lett.
88
,
207402
(
2002
).
29.
K. J.
Kubarych
,
C. J.
Milne
,
S.
Lin
,
V.
Astinov
, and
R. J. D.
Miller
,
J. Chem. Phys.
116
,
2016
(
2002
).
30.
C. J.
Milne
,
Y. L.
Li
,
T. l. C.
Jansen
,
L.
Huang
, and
R. J. D.
Miller
,
J. Phys. Chem. B
110
,
19867
(
2006
).
31.
Y. L.
Li
,
L.
Huang
,
R. J.
Dwayne Miller
,
T.
Hasegawa
, and
Y.
Tanimura
,
J. Chem. Phys.
128
,
234507
(
2008
).
32.
H.
Frostig
,
T.
Bayer
,
N.
Dudovich
,
Y. C.
Eldar
, and
Y.
Silberberg
,
Nat. Photonics
9
,
339
(
2015
).
33.
A.
Tokmakoff
,
M. J.
Lang
,
D. S.
Larsen
,
G. R.
Fleming
,
V.
Chernyak
, and
S.
Mukamel
,
Phys. Rev. Lett.
79
,
2702
(
1997
).
34.
A.
Tokmakoff
and
G. R.
Fleming
,
J. Chem. Phys.
106
,
2569
(
1997
).
35.
A.
Ma
and
R. M.
Stratt
,
J. Chem. Phys.
116
,
4962
(
2002
).
36.
A.
Ma
and
R. M.
Stratt
,
J. Chem. Phys.
116
,
4972
(
2002
).
37.
S.
Saito
and
I.
Ohmine
,
J. Chem. Phys.
125
,
084506
(
2006
).
38.
T.
Steffen
,
J. T.
Fourkas
, and
K.
Duppen
,
J. Chem. Phys.
105
,
7364
(
1996
).
39.
J. T.
Fourkas
, in
Ultrafast Infrared and Raman Spectroscopy
, edited by
M. D.
Fayer
(
Marcel Dekker
,
New York
,
2001
).
40.
X.
Sun
,
J. Chem. Phys.
151
,
194507
(
2019
).
41.
S.
Mukamel
,
V.
Khidekel
, and
V.
Chernyak
,
Phys. Rev. E
53
,
R1
(
1996
).
42.
M.
Kryvohuz
and
J.
Cao
,
Phys. Rev. Lett.
95
,
180405
(
2005
).
43.
M.
Kryvohuz
and
J.
Cao
,
Phys. Rev. Lett.
96
,
030403
(
2006
).
44.
R. M.
Stratt
,
Acc. Chem. Res.
28
,
201
(
1995
).
45.
A.
Ma
and
R. M.
Stratt
,
J. Chem. Phys.
119
,
8500
(
2003
).
46.
A.
Ma
and
R. M.
Stratt
,
Bull. Korean Chem. Soc.
24
,
1126
(
2003
).
47.
R. A.
Denny
and
D.
Reichman
,
Phys. Rev. E
63
,
065101
(
2001
).
48.
J.-Y.
Jo
,
H.
Ito
, and
Y.
Tanimura
,
Chem. Phys.
481
,
245
(
2016
).
49.
H.
Ito
,
J.-Y.
Jo
, and
Y.
Tanimura
,
Struct. Dyn.
2
,
054102
(
2015
).
50.
T.
Hasegawa
and
Y.
Tanimura
,
J. Chem. Phys.
125
,
074512
(
2006
).
51.
R.
DeVane
,
C.
Ridley
,
B.
Space
, and
T.
Keyes
,
J. Chem. Phys.
119
,
6073
(
2003
).
52.
R.
DeVane
,
C.
Ridley
,
B.
Space
, and
T.
Keyes
,
Phys. Rev. E
70
,
050101
(
2004
).
53.
R.
DeVane
,
C.
Ridley
,
B.
Space
, and
T.
Keyes
,
J. Chem. Phys.
123
,
194507
(
2005
).
54.
R.
DeVane
,
C.
Kasprzyk
,
B.
Space
, and
T.
Keyes
,
J. Phys. Chem. B
110
,
3773
(
2006
).
55.
T. E.
Markland
and
M.
Ceriotti
,
Nat. Rev. Chem.
2
,
0109
(
2018
).
56.
D.
Chandler
and
P. G.
Wolynes
,
J. Chem. Phys.
74
,
4078
(
1981
).
57.
D. M.
Ceperley
,
Rev. Mod. Phys.
67
,
279
(
1995
).
58.
M. E.
Tuckerman
,
Statistical Mechanics: Theory and Molecular Simulation
(
Oxford University Press
,
New York
,
2010
).
59.
J.
Cao
and
G. A.
Voth
,
J. Chem. Phys.
101
,
6168
(
1994
).
60.
S.
Jang
and
G. A.
Voth
,
J. Chem. Phys.
111
,
2371
(
1999
).
61.
E.
Geva
,
Q.
Shi
, and
G. A.
Voth
,
J. Chem. Phys.
115
,
9209
(
2001
).
62.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
121
,
3368
(
2004
).
63.
I. R.
Craig
and
D. E.
Manolopoulos
,
J. Chem. Phys.
122
,
084106
(
2005
).
64.
T. F.
Miller
and
D. E.
Manolopoulos
,
J. Chem. Phys.
123
,
154504
(
2005
).
65.
M.
Ceriotti
,
M.
Parrinello
,
T. E.
Markland
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
133
,
124104
(
2010
).
66.
M.
Rossi
,
M.
Ceriotti
, and
D. E.
Manolopoulos
,
J. Chem. Phys.
140
,
234116
(
2014
).
67.
B.
Leimkuhler
,
D. T.
Margul
, and
M. E.
Tuckerman
,
Mol. Phys.
111
,
3579
(
2013
).
68.
R.
Korol
,
N.
Bou-Rabee
, and
T. F.
Miller
 III
,
J. Chem. Phys.
151
,
124103
(
2019
).
69.
K. A.
Jung
,
P. E.
Videla
, and
V. S.
Batista
,
J. Chem. Phys.
148
,
244105
(
2018
).
70.
P. A. M.
Dirac
,
The Principles of Quantum Mechanics
, 4th ed. (
Clarendon Press
,
Oxford
,
1958
), Chap. 4, pp.
84
89
.
71.
D. A.
McQuarrie
,
Statistical Mechanics
(
University Science Books
,
Sausalito, CA
,
2000
), Chap. 21, pp.
507
509
.
72.
T. l. C.
Jansen
,
K.
Duppen
, and
J. G.
Snijders
,
Phys. Rev. B
67
,
134206
(
2003
).
73.
D. R.
Reichman
,
P.-N.
Roy
,
S.
Jang
, and
G. A.
Voth
,
J. Chem. Phys.
113
,
919
(
2000
).
74.
T. J. H.
Hele
,
M. J.
Willatt
,
A.
Muolo
, and
S. C.
Althorpe
,
J. Chem. Phys.
142
,
134103
(
2015
).
75.
K. A.
Jung
,
P. E.
Videla
, and
V. S.
Batista
,
J. Chem. Phys.
151
,
034108
(
2019
).
76.
K. A.
Jung
,
P. E.
Videla
, and
V. S.
Batista
, “
Ring-polymer, centroid and mean-field approximations to multi-time Matsubara dynamics
” (unpublished).
77.
J. E.
Jones
and
S.
Chapman
,
Proc. R. Soc. London Ser. A.
106
,
463
(
1924
).
78.
R. G.
Della Valle
and
E.
Venuti
,
Phys. Rev. B
58
,
206
(
1998
).
79.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Clarendon Press
,
Oxford
,
1987
), Chap. 3.
80.
X.
Sun
, Ph.D. thesis,
Brown University
,
2014
.
81.
L. C.
Geiger
and
B. M.
Ladanyi
,
J. Chem. Phys.
89
,
6588
(
1988
).
82.
J. E.
Rice
,
P. R.
Taylor
,
T. J.
Lee
, and
J.
Almlöf
,
J. Chem. Phys.
94
,
4972
(
1991
).
83.
P. A.
Madden
, in
Liquids, Freezing, and the Glass Transition
, edited by
J. P.
Hansen
,
D.
Levesques
, and
J.
Zinn-Justin
(
North-Holland
,
Amsterdam
,
1991
), pp.
547
627
.
84.
X.
Sun
and
R. M.
Stratt
,
Phys. Chem. Chem. Phys.
14
,
6320
(
2012
).
85.
X.
Sun
and
R. M.
Stratt
,
J. Chem. Phys.
139
,
044506
(
2013
).
86.
R. L.
Murry
and
J. T.
Fourkas
,
J. Chem. Phys.
107
,
9726
(
1997
).
87.
J.
Cao
,
S.
Yang
, and
J.
Wu
,
J. Chem. Phys.
116
,
3760
(
2002
).
88.
M. F.
Herman
,
E. J.
Bruskin
, and
B. J.
Berne
,
J. Chem. Phys.
76
,
5150
(
1982
).
89.
T. F.
Miller
 III
and
D. E.
Manolopoulos
,
J. Chem. Phys.
122
,
184503
(
2005
).
90.
T. E.
Markland
,
J. A.
Morrone
,
B. J.
Berne
,
K.
Miyazaki
,
E.
Rabani
, and
D. R.
Reichman
,
Nat. Phys.
7
,
134
(
2011
).
91.
T. E.
Markland
,
J. A.
Morrone
,
K.
Miyazaki
,
B. J.
Berne
,
D. R.
Reichman
, and
E.
Rabani
,
J. Chem. Phys.
136
,
074511
(
2012
).
92.
G. R.
Medders
,
A. W.
Götz
,
M. A.
Morales
,
P.
Bajaj
, and
F.
Paesani
,
J. Chem. Phys.
143
,
104102
(
2015
).
93.
X.
Sun
,
B. M.
Ladanyi
, and
R. M.
Stratt
,
J. Chem. Phys. B
119
,
9129
(
2015
).
94.
X.
Sun
,
P.
Zhang
,
Y.
Lai
,
K. L.
Williams
,
M. S.
Cheung
,
B. D.
Dunietz
, and
E.
Geva
,
J. Phys. Chem. C
122
,
11288
(
2018
).
95.
J.
Savolainen
,
S.
Ahmed
, and
P.
Hamm
,
Proc. Natl. Acad. Sci. U. S. A.
110
,
20402
(
2013
).
96.
D.
Sidler
and
P.
Hamm
,
J. Chem. Phys.
150
,
044202
(
2019
).
97.
A.
Berger
,
G.
Ciardi
,
D.
Sidler
,
P.
Hamm
, and
A.
Shalit
,
Proc. Natl. Acad. Sci. U. S. A.
116
,
2458
(
2019
).
98.
J. R.
Cendagorta
,
Z.
Bačić
, and
M. E.
Tuckerman
,
J. Chem. Phys.
148
,
102340
(
2018
).
99.
X.
Sun
and
E.
Geva
,
J. Chem. Theory Comput.
12
,
2926
(
2016
).
100.
X.
Sun
and
E.
Geva
,
J. Chem. Phys.
145
,
064109
(
2016
).
101.
W.
Kuehn
,
K.
Reimann
,
M.
Woerner
,
T.
Elsaesser
,
R.
Hey
, and
U.
Schade
,
Phys. Rev. Lett.
107
,
067401
(
2011
).

Supplementary Material