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.

## I. INTRODUCTION

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 dynamics^{1–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 theorem^{8} 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 3*N* 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 spectroscopy^{15–23} introduces an interaction (*H*_{int} = −**E**_{1} ⋅ **Π** ⋅ **E**_{2}/2, where **E**_{1/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 CS_{2},^{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}

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/*k*_{B}*T* (*k*_{B} is Boltzmann’s constant) can be considered as a propagator *e*^{−iĤ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 3*N*-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.

## II. THEORY

### A. 2D Raman spectroscopy: Classical molecular dynamics approaches

The optical observable in fifth-order 2D Raman spectroscopy can be expressed in terms of the two-time response function as^{9}

Here, $\Pi ^$ represents the many-body polarizability tensor, *Â*(*t*) = *e*^{iĤt/ℏ}*Âe*^{−iĤt/ℏ} is the operator evolved in the Heisenberg picture, $[\xc2,B^]=\xc2B^\u2212B^\xc2$ is the quantum commutator, and $\rho ^eq=e\u2212\beta \u0124/Z$ is the equilibrium density operator with system Hamiltonian *Ĥ*, inverse temperature *β* = 1/*k*_{B}*T*, and partition function *Z* = Tr[*e*^{−βĤ}]. By selecting the polarization of the perturbing electric fields **E**_{i} 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 $\Pi ^=\Pi ^=\Pi ^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}

where

Here, *A* = *A*(**r**, **p**) and *B* = *B*(**r**, **p**) represent classical dynamical variables, **r** = (**r**_{1}, …, **r**_{N}) are the 3*N*-dimensional coordinates, and **p** = (**p**_{1}, …, **p**_{N}) are the corresponding momenta of an *N*-atom system. Using this quantum–classical correspondence and noting that {*A*, *ρ*_{eq}} = −*β*$\u0226$*ρ*_{eq} with $\u0226$ = d*A*/d*t*,^{71} the 2D Raman response function can then be expressed as^{35}

where ⟨·⟩ = ∫ d**r**d**p***ρ*_{eq} represents a classical ensemble average with Hamiltonian *H*, phase-space density *ρ*_{eq} = *e*^{−βH}/*Z*, and classical partition function *Z* = ∫ d**r**d**p***e*^{−β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., $\Pi ^=\Pi ^(r)$. The Poisson bracket in the classical 2D response function can be expressed as

and the classical 2D response is thus

Here, **J**^{rp}(*t*) = *∂***r**_{t}/*∂***p**_{0} is the stability matrix (or Jacobian matrix) satisfying the equations of motion $J\u0308rp(t)=\u2212(1/m)D(t)\u22c5Jrp(t)$ with initial condition **J**^{rp}(0) = **0**, mass *m*, and dynamical matrix (or instantaneous Hessian) $D(t)=\u2207\u2207V(r)rt$.^{35} The stability matrix describes the sensitivity of positions at time *t*, **r**_{t}, if the initial momenta, **p**_{0}, 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,

Here, Π_{±Π(0)} (*t*_{2}) represents the many-body polarizability at time *t*_{2} calculated from nonequilibrium trajectories that are generated by an ultrafast Raman perturbation *H*_{int} = ±(−**E**_{1} ⋅ **Π** ⋅ **E**_{2}/2) at time *t* = 0, with external electric fields **E**_{1/2} of time duration Δ*t*.^{40,72} The extra forces due to the Raman perturbation are given by Δ**F** = −$\u2207$_{r}*H*_{int} = ±**E**_{1} ⋅ $\u2207$_{r}**Π** ⋅ **E**_{2}/2. The 2D Raman response function in the hybrid equilibrium–nonequilibrium is then given by

where $\Pi \u0307(\u2212t1)$ 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

where *C*(*t*_{1}, *t*_{2}) represents the classical two-time correlation function of the fluctuations of the many-body polarizability, namely,

with $\delta \Pi =\Pi \u2212\Pi $. 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 function^{54} 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}

### B. 2D Raman spectroscopy: Ring-polymer molecular dynamics approach

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

where the *Q* factors are defined as

with

and $K\u0303sym(\omega 1,\omega 2)$ and $K\u0303asym(\omega 1,\omega 2)$ are the double Fourier transforms of the symmetric (real) and asymmetric (imaginary) DKT correlation functions, namely,

Here, the DKT correlation function is defined as^{69,73}

It is noted that in Ref. 69, we used a different choice of time variables *t*_{1}, *t*_{2} in the response function. However, with a change of variables (*t*_{1}, *t*_{2}) → (*t*_{1}, *t*_{1} + *t*_{2}) 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, as^{69}

Moreover, semiclassical methodologies including multi-time Matsubara dynamics^{74,75} and multi-time RPMD^{62,69} were developed to evaluate the symmetric DKT correlation function *K*^{sym}(*t*_{1}, *t*_{2}).^{69,75} In this study, we use RPMD to approximate the symmetric DKT in terms of polarizability fluctuations as

where

Here, *f* = 3*NP* 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(\alpha )})(j=1,\u2026,N;\alpha =1,\u2026,P)$ are the ring-polymer coordinates (with a similar definition for the conjugate momenta **P**). *Z*_{P} is the *P*-bead path-integral representation of the partition function,

*A*_{P} is the ring-polymer representation of any observable *A*,

*H*_{P} is the ring-polymer Hamiltonian,

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

with *ω*_{P} = 1/(*β*_{P}*ℏ*) and $rj(0)\u2261rj(P)$. In Eq. (19), **R**_{t} denotes the ring-polymer coordinates at time *t* evolved from **R**_{0} at time *t* = 0 using the Hamiltonian *H*_{P}.

The RPMD approximation [Eq. (18)] can be shown to be an exact path-integral discretization of the symmetric DKT in the limit *t*_{1}, *t*_{2} → 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}

## III. COMPUTATIONAL DETAILS

### A. Lennard-Jones model

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

where *r*_{ij} = |**r**_{ij}| = |**r**_{i} − **r**_{j}| 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 *ε*/*k*_{B} = 222 K, *σ* = 4.099 Å, and mass *m* = 131.293 amu (with reduced time unit $\tau LJ=m\sigma 2/\epsilon =3.47$ ps) at the phase point *T*^{*} ≡ *Tk*_{B}/*ε* = 1.00, *ρ*^{*} ≡ *ρσ*^{3} = 0.80 and (ii) a LJ fluid representing neon with parameters *ε*/*k*_{B} = 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, $\u210f*\u2261\u210f/m\u03f5\sigma 2$) that is dependent on the mass of the particle (*ℏ*^{*} ≈ 0.009 95 for Xe and *ℏ*^{*} ≈ 0.094 53 for Ne).^{78}

### B. Molecular dynamics simulation details

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 integrator^{79} 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 × 10^{8} liquid configurations from 500 independent trajectories of 1.73 ns were averaged. In the hybrid equilibrium–nonequilibrium MD simulation, 1 × 10^{8} liquid configurations generated from equilibrium MD were sampled to generate 1 × 10^{8} pairs of positively and negatively perturbed nonequilibrium trajectories by applying extra forces with external electric fields *E*_{1} = *E*_{2} = 5.0 V/Å. In the DRSK TCF calculation, 1 × 10^{8} liquid configurations from 500 independent trajectories of 1.73 ns were averaged to obtain $\delta \Pi (0)\delta \Pi t1\delta \Pi 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 × 10^{8} 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 |*t*_{1,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 × 10^{8} ring-polymer configurations from 500 independent trajectories (2.27 ns each) were averaged to obtain $\delta \Pi (0)\delta \Pi t1\delta \Pi t1+t2P$. The TCF was computed with a maximum *t*_{1} and *t*_{2} of 2.04 ps. Table I summarizes the simulation parameters.

. | N
. | P
. | δt/τ_{LJ}
. | Traj. . | Config. . | β⟨KE⟩/N
. |
---|---|---|---|---|---|---|

Xe (RPMD) | 108 | 16 | 0.0005 | 200 | 4 × 10^{8} | 1.50 ± 0.03 |

Xe (MD) | 108 | 1 | 0.0005 | 500 | 1 × 10^{8} | 1.50 ± 0.07 |

Ne (RPMD) | 64 | 64 | 0.001 | 500 | 1 × 10^{8} | 1.85 ± 0.04 |

Ne (MD) | 108 | 1 | 0.001 | 500 | 1 × 10^{8} | 1.50 ± 0.08 |

. | N
. | P
. | δt/τ_{LJ}
. | Traj. . | Config. . | β⟨KE⟩/N
. |
---|---|---|---|---|---|---|

Xe (RPMD) | 108 | 16 | 0.0005 | 200 | 4 × 10^{8} | 1.50 ± 0.03 |

Xe (MD) | 108 | 1 | 0.0005 | 500 | 1 × 10^{8} | 1.50 ± 0.07 |

Ne (RPMD) | 64 | 64 | 0.001 | 500 | 1 × 10^{8} | 1.85 ± 0.04 |

Ne (MD) | 108 | 1 | 0.001 | 500 | 1 × 10^{8} | 1.50 ± 0.08 |

### C. Spectroscopic observable: Many-body polarizability

The many-body polarizability **Π** needed for the calculation of the fifth-order full-polarized 2D Raman response *R*^{(5)} (*t*_{1}, *t*_{2}) 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*),

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

with

where $r^\u2032=r\u2032/r\u2032,\u2003r\u2032=|r\u2032|,and\u2003(rjk=rk\u2212rj)$.

In the case of LJ atomic liquids, each atomic polarizability tensor *α*_{j} = *α*_{j}**1** 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

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}

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

In our case, *A*_{zz}, *B*_{zz}, and *C*_{zz} 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.

## IV. RESULTS

### A. Nuclear quantum effects in Lennard-Jones liquids

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 by^{88}

where $r\xafj=1P\u2211\alpha =1Prj(\alpha )$ 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\xafj$, 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, ⟨*KE*⟩_{1} = 3*N*/(2*β*). Thus, values of *β*⟨*KE*⟩_{P}/*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.

As a final descriptor of NQE on LJ liquids, we performed RPMD simulations to compute the velocity autocorrelation function defined as^{64,89}

which represents an approximation to the Kubo transformed velocity autocorrelation function. In the previous equation, $v\xafj=1P\u2211\alpha =1Pvj(\alpha )$ 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 *C*_{vv}(*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 ^{4}He droplets below 1 K.^{57}

### B. 2D Raman spectroscopy of liquid Xe

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 $\delta \Pi (0)\delta \Pi t1\delta \Pi 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*(*t*_{1}, *t*_{2}) = *f*(−*t*_{1}, −*t*_{2})] ^{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).

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 *t*_{1} ∼ 50 fs, *t*_{2} ∼ 340 fs, in consistency with previous results.^{35,36} Note also the lack of an echo signal along the diagonal *t*_{1} = *t*_{2}^{35,36,45,46} and that *R*^{(5)} (*t*_{1}, *t*_{2}) = 0 along the *t*_{2} = 0 axis.

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 *t*_{1} ∼ 80 fs, *t*_{2} ∼ 340 fs, and a response that vanishes (within numerical error) along the *t*_{2} = 0 axis. However, both DKT and DRSK responses show a faster decay along the *t*_{2} axis for *t*_{2} < 0.4 ps, and a different spectral decay in the region *t*_{1} < 0.5 ps, *t*_{2} < 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}

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 *t*_{2} = 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 *t*_{2} < 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.

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, *K*^{sym}(*t*_{1}, *t*_{2}) → *C*(*t*_{1}, *t*_{2}).^{69} Moreover, it is straightforward to show that in this high temperature limit (see the Appendix),

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

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

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.

### C. 2D Raman spectroscopy of liquid Ne

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 $\delta \Pi (0)\delta \Pi t1\delta \Pi 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 *t*_{1} = 50 fs, *t*_{2} = 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 *t*_{2} 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).

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 *t*_{2} = 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.

## V. CONCLUDING REMARKS

In the present work, we have shown that the recently developed double Kubo transform (DKT) approach^{69} 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*-H_{2},^{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 methods^{98–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 THz^{101} 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.

## SUPPLEMENTARY MATERIAL

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

## DATA AVAILABILITY

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

## ACKNOWLEDGMENTS

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.

### APPENDIX: HIGH TEMPERATURE LIMIT OF $Q\xb1$

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 *Q*_{1} [Eq. (13)],

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