# Analysis of 2D THz-Raman spectroscopy using a non-Markovian Brownian oscillator model with nonlinear system-bath interactions

###### Abstract

We explore and describe the roles of inter-molecular vibrations employing a Brownian oscillator (BO) model with linear-linear (LL) and square-linear (SL) system-bath interactions, which we use to analyze two-dimensional (2D) THz-Raman spectra obtained by means of molecular dynamics (MD) simulations. In addition to linear absorption (1D IR), we calculated 2D Raman-THz-THz, THz-Raman-THz, and THz-THz-Raman signals for liquid formamide, water, and methanol using an equilibrium non-equilibrium hybrid MD simulation. The calculated 1D IR and 2D THz-Raman signals are compared with results obtained from the LL+SL BO model applied through use of hierarchal Fokker-Planck equations with non-perturbative and non-Markovian noise. We find that all of the qualitative features of the 2D profiles of the signals obtained from the MD simulations are reproduced with the LL+SL BO model, indicating that this model captures the essential features of the inter-molecular motion. We analyze the fitted 2D profiles in terms of anharmonicity, nonlinear polarizability, and dephasing time. The origins of the echo peaks of the librational motion and the elongated peaks parallel to the probe direction are elucidated using optical Liouville paths.

###### pacs:

Valid PACS appear here^{†}

^{†}preprint: APS/123-QED

## I Introduction

Molecular vibrations in condensed phases play an essential role in various dynamic processes, including inter- and intra-molecular couplings, and solvent dynamics, all of which entail energy exchange as well as thermal excitations and relaxations.Mukamel-Book () Multidimensional vibrational spectroscopy techniques make it possible to experimentally distinguish such processes due to the sensitivity of the nonlinear response functions utilized in these techniques to complex dynamics.Cundiff-PT-2013-7 (); Mukamel-ACR-2009-42 () For intra-molecular vibrations, the roles of relaxation and dephasing are well understood both theoretically and experimentally due to the advent of infrared (IR) laser technologies. Methods of analysis with theoretical models that utilize molecular dynamic simulations have also been developed to elucidate multidimensional IR signals.Asbury-JPCA-2004-108 (); Auer-PNAS-2007-104 () Because the primary inter-molecular modes, which are the objects of study in 2D IR spectroscopy, can be separated from the other modes, as in the case of the OH stretching mode in liquid water, stochastic models whose parameters are obtained from classical molecular dynamics (MD) simulations have been useful for analysis of the inter-molecular vibrational modes. For inter-molecular vibrational modes, two-dimensional Raman spectroscopyTanimura-JCP-1993-99 () was for a long time the only two-dimensional spectroscopy that could be used for experimental study. However, due to technical difficulties, such investigations have been carried out only for ,Kaufman-PRL-2002-88 (); Kubarych-IRPC-2003-22 (); Kubarych-JCP-2002-116 (); Kubarych-CPL-2003-369 () Benzene,Milne-JPCB-2006-40 () and formamideLi-JCP-2008-128 () liquids. Theoretical investigations have also been limited, due to the availability of experimental data and limitations on computational power for simulations.

Two-dimensional THz-Raman spectroscopy, which has been studied both theoreticallyHamm-JCP-2012-136 (); Hamm-JCP-2012-136-note (); Hamm-JCP-2014-141 (); Ito-JCP-2014-141 () and experimentally,Savolainen-PNAS-2013-110 () has created a new possibility for investigating the details of inter-molecular vibrations. In 2D Raman spectroscopy, the observable is defined in terms of the three-body response function for the polarizability of the system, as , where represents the thermal average and is the Heisenberg operator for an arbitrary operator .Tanimura-JCP-1993-99 () In the case of 2D THz-Raman spectroscopy, the response function consists of one polarizability, , and two dipole operators, , and there are three different measurements, which depend upon the sequence of the Raman and THz pulses as , , and . While inter-molecular vibrational modes are usually both Raman and IR active, the types of information that we can obtain from the 2D Raman signal and each of three THz-Raman signals are different, due to the role of the nonlinear polarizability. Because each of the above-mentioned response functions is defined in terms of the three-body correlation function, the signal will vanish if the system is harmonic and if the total dipole moment and polarizability are linear functions of the collective coordinate, , representing the inter-molecular vibration, because there is an odd number of Gaussian integrals involved in the response function: . The dipole moment is approximated reasonably well as a linear function of as , because the total dipole moment is a linear function of the distance between the charges in the system, and the nonlinear dipole-induced dipole interactions are weak. However, the contribution of the non-linear polarization is not negligible, because the polarizability originates in the electronic states of molecules, which depend on the complex configurations of the atoms and molecules. For this reason the polarizability is expressed in a Taylor expansion form as . Because has this non-linear form, the three response functions give above, representing the observables in 2D THz-Raman spectroscopy experiments, provide information about three different physical processes. Hamm-JCP-2012-136-note () Contrastingly, because 2D Raman spectroscopy experiments measure just a single observables, they do not provide such a detailed picture of the physical system. The richness of the information obtained through 2D THz-Raman spectroscopy allows for a detailed analysis of inter-molecular vibrational modes. Note that the optical setup for the TTR measurement differs significantly from those for the RTT and TRT measurements. This is because the RTT and TRT responses are detected as the emission of THz signals, while the TTR response is detected as an induced Raman signal.

Although we can obtain relatively reliable 2D THz-Raman signals using the full MD simulation techniques developed for 2D Raman spectroscopy,Ma-PRL-2000-85 (); Saito-PRL-2002-88 (); Saito-JCP-2003-119 (); Saito-JCP-2006-125 (); Nagata-JCP-2007-126 (); Jansen-JCP-2001-114 (); Jansen-JCP-2003-63 (); Hasegawa-JCP-2006-125 () analysis of the spectra is not straightforward, due to the complexity of the 2D profiles of the signals, which arises from the complexity of the inter-molecular vibrational modes. As demonstrated by 2D Raman and 2D IR spectroscopy studies, a model-based analysis is useful for treating this problem, because the 2D profile of the signal is so sensitive to the underlying dynamics that the complex 2D profile cannot be reproduced without capturing the essential features of the vibrational modes.Tanimura-ACR-2009-42 () While stochastic models, which can be regarded as Brownian models with non-linear system-bath interactions,Ishizaki-JCP-2006-125 (); Ishizaki-JPC-2007-111 () are recognized as versatile models for analyzing intra-molecular modes observed in 2D IR spectroscopy experiments, it is not clear if such models are useful in the 2D THz-Raman case. This is because in contrast to the intra-molecular modes, which are clearly definable with the normal mode picture, the inter-molecular modes are not localized and change in time due to changes in the configuration of the system molecules.

In this paper, we explore the possibility of characterizing inter-molecular modes using a Brownian model with linear-linear (LL) and square-liner (SL) interactions utilizing 2D THz-Raman signals obtained from MD simulations. In order to treat a non-perturbative, non-Markovian, and nonlinear system-bath interaction, which is necessary to describe the effects of homogeneous and inhomogeneous broadening in a unified manner, we employ the hierarchal equations of motion approach.Tanimura-JPSJ-1989-58 (); Tanimura-PRA-1991-43 (); Tanimura-JCP-1992-96 (); Steffen-JPSJ-2000-69 (); Tanimura-JPSJ-2000-69 (); Kato-JCP-2002-117 (); Kato-JCP-2004-120 (); Tanimura-JPSJ-2006-75 (); Tanimura-JCP-accepted () The properties of inter-molecular motion are investigated using the fitted model.

This paper is organized as follows. In Sec. II, we explain the methodology for calculating 2D THz-Raman signals from full MD simulations. In Sec. III, we present the LL+SL BO model and the hierarchal equations of motion formalism. We then show how this formalism can be used to calculate 2D signals. The MD and fitted results obtained from the LL+SL BO model are presented and analyzed in Sec. IV. Section V is devoted to concluding remarks.

## Ii Full MD simulation

While, to this time, the experimentally obtained 2D THz-Raman signals are limited to the case of liquid water, in this paper we analyze 2D signals obtained from full MD simulations for formamide, water, and methanol. We chose these liquids from among many substances that have been investigated in full MD studies of the 2D Raman and 2D THz-Raman spectroscopy as characteristic examples of 2D THz-Raman signals. The MD simulation results used in the present study of these molecules for 2D Raman-THz-THz (RTT) and THz-Raman-THz (TRT) signals were originally presented in a previous study.Ito-JCP-2014-141 () Nevertheless, here we repeated the full MD simulations in order to also obtain THz-THz-Raman (TTR) and infrared absorption signals, in addition to the RTT and TRT signals. Moreover, we employed the Ewald sum for the evaluation of the dipole and polarizability, in addition to the force fields. This contrasts with the situation in previous studies, in which only force fields were computed with the Ewald sum. The change in the resulting signals due to the use of the Ewald sum in the computation of the dipole and polarizability, however, are small.

### ii.1 Models and simulation details

Based on the MD simulations, we calculated the linear absorption (1DIR) spectrum and 2D THz-Raman signals of liquid formamide, water, and methanol. Each system consisted of 108 molecules in a cubic box with periodic boundary conditions. The interactions between the molecules were modeled by a modified T potential,Sagarik-JCP-1987-86 (); Wojcik-JCP-2000-113 () the TIP4P/2005 potential,Abascal-JCP-2005-23 () and the B3 potentialWalser-JCP-2000-112 () for formamide, water, and methanol, respectively.

The interaction potentials were cut off smoothly at a distance equal to a half the length of the system using a switching function, and the long-range Coulomb interactions were calculated with the Ewald sum. The intra-molecular geometries were kept rigid throughout the simulations, using a constraint provided by the RATTLE algorithm. The equations of motion were integrated using the velocity-Verlet algorithm with time steps of for formamide and for water and methanol. The system volume and total energy were fixed after the completion of the isothermal simulations carried out for equilibration. The conditions of the simulation were set such that the average densities were for formamide, for water, and for methanol. The temperature was set to . The permanent molecular polarizability of each liquid was utilized with the atomic polarizability for formamide and methanol,Applequist-JACS-1972-94 () and the Huiszoon polarizability for water.Huiszoon-MP-1986-58 ()

### ii.2 Molecular polarizability and dipole

While the MD simulations were carried out using the permanent polarizability, we calculated the 2D THz-Raman signals using a full-order dipole-induced-dipole (DID) polarizability model. We did this because the 2D profiles are extremely sensitive to the accuracy of the calculated optical observables. In the DID polarizability model, the expression determining the polarizability of a molecule includes contributions from other molecules, and interactions between molecules are defined with respect to the centers of individual molecules. Ladanyi-CPL-1985-121 (); Applequist-JACS-1972-94 () The total polarizability of the system in a MD simulation is given by , where is the polarizability of the th molecule, expressed as

(1) |

In this expression, is the permanent molecular polarizability of the th molecule in isolation in the laboratory frame, and is the dipole-dipole interaction tensor

(2) |

Here, is the vector from the center of mass of molecule to the center of mass of molecule , and . Also, and are the unit matrix and the tensor product, respectively. In order to properly take into account the effect of the long-range interaction on the molecular polarizability, we employ the Ewald sum for . This effect has been ignored in previous MD simulations of 2D Raman and 2D THz-Raman spectroscopy systems. However, the contribution of this effect in the 2D THz-Raman case is minor in comparison with that in the 2D Raman case.

The total dipole moment is evaluated as , where and are the permanent and induced molecular dipole moments of molecule , respectively. The induced molecular dipole moment is expressed in terms of the interaction tensor as

(3) |

where is the electrostatic field at molecule created by all the other molecules in the system. This is evaluated as , where is the vector between the center of mass of molecule and that of atom in molecule .

### ii.3 One- and two-dimensional signals

It should be noted that the majority of MD simulations of 2D IR spectroscopy systems performed to this time have been carried out to obtain the parameter values for stochastic models. Full MD simulations have been carried out mostly in the cases of low frequency vibrational modes.Hasegawa-JCP-2008-128 (); Yagasaki-JCP-2008-128 (); Jeon-JPB-2014-118 () This is because the primary inter-molecular modes, which are the objects of study in 2D IR spectroscopy, can be separated from the other modes rather easily, as in the case of the OH stretching mode in liquid water. Contrastingly in the 2D THz-Raman case, it is not easy to find primary modes, because the objects of study in this cases are inter-molecular vibrations that depend on the complicated nature of molecular ensembles, whose configurations change in time. Thus, we have to evaluate 2D signals directly from the MD simulations. Because quantum mechanical effects are minor for low-frequency inter-molecular modes, due to their small thermal activation energies, unlike the case of intra-molecular motion, Sakurai-JPCA-2011-115 () and because 2D THz-Raman spectroscopy employs the three-body correlation function with two time variables, instead of the four-body correlation function with three time variables employed in 2D IR spectroscopy, the full MD simulation approach is practical. For this reason with our approach, we were able to carry out full MD simulations to directly evaluate 2D signals.

Although our MD and model calculations are fully classical, we start from the quantum expressions for the response functions, because their classical expressions in the MD and model calculations are most easily derived by taking the classical limit of the quantum expressions. The optical observables in 1D and 2D spectroscopies are represented respectively expressed by two- and three-body response functions of the formsTanimura-JCP-1993-99 (); Tanimura-JPSJ-2006-75 ()

(4) |

and

(5) |

where , , and can be the total dipole moment, , or the total polarizability of the molecules, . For low-frequency inter-molecular vibrations, we can take the classical limit, . The commutator and operators are then replaced by the Poisson bracket and c-number observables as

(6) |

Using for a molecular Hamiltonian , we obtain the expression for the¡ linear response function, for example, for 1D IR as

(7) |

and

(8) |

It should be noted that the quantum correction factor, , is usually included in Eq. (8) to allow comparison with experimentally obtained IR signals, but here we do not include it. Instead, as the classical limit, we only multiply by . Berens-JCP-1981-74 () We can evaluate the above quantities easily by calculating from samples of molecular trajectories that are obtained from the equilibrium MD simulation.

For the 2D case, the response function in the classical limit is expressed asMa-PRL-2000-85 (); Saito-PRL-2002-88 ()

(9) |

As in the 1D case, we can calculate the above response function using and evaluated from the molecular trajectories and obtained from the equilibrium MD simulations. Ma-PRL-2000-85 (); Saito-PRL-2002-88 (); Saito-JCP-2003-119 (); Saito-JCP-2006-125 (); Nagata-JCP-2007-126 () The convergence of the signal is, however, very slow due to the effect of the stability matrix element in the double Poisson brackets. However, there is different approach, the non-equilibrium finite field approach, that does not have this convergence problem. In this approach, the double Poisson brackets are evaluated as , where is the observable corresponding to calculated from the trajectories subjected to the weak perturbations , with the electric field acting on .Jansen-JCP-2001-114 (); Jansen-JCP-2003-63 () But, this approach is computationally intensive, and for this reason, here we employed a hybrid approach, which utilizes both the equilibrium and non-equilibrium approaches in order to reduce the computational cost further.Hasegawa-JCP-2006-125 ()

In our hybrid approach, we evaluate with equilibrium MD simulations, while with non-equilibrium MD simulations. As a result, the hybrid expressions for the 2D Raman-THz-THz, THz-Raman-THz, and THz-THz-Raman signals become

(10) | ||||

(11) |

and

(12) |

where and are the external electric field of the th pulse and the inverse temperature divided by the Boltzmann constant.

## Iii Model calculation

### iii.1 Brownian oscillator model with nonlinear interaction

In order to analyze the 2D signals of inter-vibrational modes, we consider a model that consists of a primary oscillator mode nonlinearly coupled to the other modes, which are regarded as a bath system. This bath system is represented by an ensemble of harmonic oscillators. The primary mode may change in time or be inhomogeneously distributed. We can describe both situations within a unified framework by adjusting the bath parameter variables. The model is constructed by extending a Brownian (or Caldeira-Leggett) HamiltonianCaldeira-PA-1983-121 (); Grabert-PR-1988-115 () to include a nonlinear system-bath interaction. We write

(13) |

where

(14) |

is the Hamiltonian for the system with mass , momentum and potential ,

(15) |

is the bath Hamiltonian with the momentum, coordinate, mass, and frequency of the th bath oscillator given by , , and , respectively, and

(16) |

is the system-bath interaction, which consists of linear-linear (LL) and square-linear (SL) system-bath interactions, , with coupling strengths , , and .Okumura-PRE-1997-56 () This model has been used to derive predictions for 2D RamanSteffen-JPSJ-2000-69 (); Tanimura-JPSJ-2000-69 (); Kato-JCP-2002-117 (); Kato-JCP-2004-120 (); Tanimura-JPSJ-2006-75 () and 2D IR signals.Ishizaki-JCP-2006-125 (); Ishizaki-JPC-2007-111 (); Tanimura-ACR-2009-42 () The last term of the bath Hamiltonian is the counter-term, which maintains the translational symmetry of the system in the case .

The sum of the bath coordinates acts as a collective coordinate that modulates the system.Tanimura-JPSJ-2006-75 () As illustrated in Ref. Tanimura-ACR-2009-42, , while the LL interaction shifts the potential, the SL interaction changes its curvature. Although in the anharmonic potential case, the LL interaction also changes the curvature of the potential, we can ignore this effect if the anharmonicity is weak.Sakurai-JPCA-2011-115 (); Tanimura-ACR-2009-42 () Next we introduce the spectral distribution function, , which characterizes the bath and system-bath coupling. We assume that has an Ohmic form with a Lorentzian cutoff:Tanimura-JPSJ-1989-58 (); Tanimura-PRA-1991-43 (); Tanimura-JCP-1992-96 (); Steffen-JPSJ-2000-69 (); Tanimura-JPSJ-2000-69 (); Kato-JCP-2002-117 (); Kato-JCP-2004-120 (); Tanimura-JPSJ-2006-75 (); Tanimura-JCP-accepted ()

(17) |

where is the system-bath coupling strength, and represents the width of the spectral distributuion.

Writing the classical collective coordinate corresponds as , we have the correlation function, . This indicates that the bath oscillators interact with the system in the form of Gaussian Markovian noise with correlation time .Tanimura-JPSJ-2006-75 ()

Because the SL interaction affects the frequency of the potential, the fast modulation limit () of the SL interaction corresponds to the case of a homogeneous distribution, as depicted in Fig. 1(a). By contrast, the slow modulation limit () corresponds to the case of an inhomogeneous distribution, as depicted in Fig. 1(b). Note that SL or LL+SL BO model has the “mode mixing in polarization”Saito-JCP-108-1998 () because the modes included in collective coordinate are frequency distributed by SL interaction and are mixed by nonlinear polarizability .

### iii.2 Classical hierarchal Fokker-Planck equations

Because we wish to explore the effects of anharmonicity, nonlinear polarizability, vibrational dephasing, and homogeneous and inhomogeneous broadening within a unified framework, we must employ a kinetic equation that can treat thermal fluctuations as well as dissipation in a non-perturbative, non-Markovian manner. The reduced hierarchal equations of motion (HEOM) satisfy all of the requirements mentioned above and is ideal for the present study.Tanimura-JPSJ-1989-58 (); Tanimura-PRA-1991-43 (); Tanimura-JCP-1992-96 (); Tanimura-JCP-accepted (); Kato-JCP-2002-117 (); Kato-JCP-2004-120 (); Tanimura-JPSJ-2006-75 (); Steffen-JPSJ-2000-69 (); Tanimura-JPSJ-2000-69 () While we must use the quantum form of the equations to calculate the signals for high frequency intra-molecular modes,Sakurai-JPCA-2011-115 () we can employ the classical form for low frequency inter-molecular modes.

For the LL+SL BO Hamiltonian, given in Eqs. (13)-(16), the HEOM for the classical distribution function are expressed asKato-JCP-2002-117 (); Kato-JCP-2004-120 (); Tanimura-JPSJ-2006-75 ()

(18) |

for and

(19) |

In the HEOM approach, only the first element, , has physical meaning, while the other elements, , are introduced in the numerical calculations in order to treat the non-perturbative, non-Markovian system-bath interaction. We choose to satisfy , where is the characteristic frequency of the system.

The classical Liouvillian of the system, , is defined by

(20) |

where the dash is defined as for an arbitrary function . The operators and describe the energy exchange between the system and the heat bath for the inverse correlation time . They are defined as

(21) |

and

(22) |

with and

(23) |

The thermal equilibrium distribution, , is expressed in terms of the HEOM elements evaluated from the steady-state solution of the HEOM. Note that Eqs. (18)-(19) reduce to the Kramers equation in the limit with .Risken-Book ()

Hereafter, we employ the dimensionless coordinate and momentum defined by and , where represents the fundamental frequency. The potential is then assumed to be

(24) |

where is the cubic anharmonicity of the potential. The other variables, , , , and , are also normalized accordingly.

### iii.3 One- and two-dimensional signals

To apply the HEOM formalism, we express the response functions in terms of the time-propagation operator. Then, Eqs. (4) and (5) can be rewritten as

(25) |

and

(26) |

where we have employed the hyperoperator defined as , is the Green’s function of the system Hamiltonian without a laser interaction, and is the equilibrium state. The above equations represent the time evolution of the system under laser excitation. For example, Eq. (26) can be interpreted as follows. The system is initially in the equilibrium state and is then modified as a result of the first laser pulse via the dipole interaction by . It then propagates for time under . The system is next excited through the second laser pulse by and propagates for time under . Finally, the expectation value of the polarizability at is generated through the laser pulses by .Tanimura-JPSJ-2006-75 ()

The classical expressions for the response functions can be obtained from the above with the use of the Wigner transformation.Tanimura-CP-1998-233 () In this case, an arbitrary operator is replaced by . For 1D IR and 2D THz-Raman spectroscopies, we have

(27) |

and

(28) | ||||

(29) |

and

(30) |

respectively. Here, is the relative intensity of the nonlinear polarizability. The Green’s function is now expressed in terms of the classical Liouvillian as , and is the equilibrium distribution. In order to apply the HEOM formalism, we express the time-dependent Wigner function , such as and , in terms of the HEOM member and , and the determine its time evolution through Eqs.(18) and (19).Steffen-JPSJ-2000-69 (); Tanimura-JPSJ-2000-69 (); Kato-JCP-2002-117 (); Kato-JCP-2004-120 (); Tanimura-JPSJ-2006-75 () In the HEOM formalism, the equilibrium distribution, , is also expressed using the HEOM elements evaluated from the steady-state solution of Eqs. (18) and (19). In the strong coupling case, we employ an eigenfunctional representation of the momentum space for numerical convenience, as discussed in Appendix A, while, in other cases, we solve Eqs. (18) and (19). Note that, as shown in Ref. Tanimura-PRA-1991-43, , the hierarchal Fokker-Planck approach is equivalent to the generalized Langevin approach. However, because the Fokker-Planck approach does not require sampling of system trajectories, unlike the Langevin approach, it is numerically advantageous, especially for calculating nonlinear response functions, for which the trajectories are unstable.

## Iv Results and Discussion

Because the 2D profiles of the signal must be constructed from complex motion in a complicated manner, the analysis of the signal profile is not straightforward. Nevertheless, properly accounting for the components of signals to these profiles allows us to perform a detailed analysis of the inter-molecular vibrational motion on the basis of both experiential and theoretical results. Model-based studies of the 2D profiles of signals are helpful to identify the underlying physical mechanisms, because it is necessary to capture the essential features of the inter-molecular motion in order to reproduce the complex 2D profile from a simple model. Analyses of this kind have employed LL + SL BO models for the 2D IRIshizaki-JCP-2006-125 (); Ishizaki-JPC-2007-111 (); Tanimura-ACR-2009-42 (); Sakurai-JPCA-2011-115 () and 2D Raman cases.Steffen-JPSJ-2000-69 (); Tanimura-JPSJ-2000-69 (); Kato-JCP-2002-117 (); Kato-JCP-2004-120 (); Tanimura-JPSJ-2006-75 () However, their applicability has not been fully explored because of the limited availability of experimental and numerical data. 2D THz-Raman measurements, which are applicable not only to the study of liquids, but also to the problem of distinguishing optical processes through use of RTT, TRT, and TTR measurements, provide the opportunity to explore the possibilities of model-based analysis of 2D profiles. Here, we examine the LL+SL BO model and use it to reproduce all three THz-Raman signal profiles obtained from full MD simulations. This is done by choosing the parameter values of the model so as to realize the best agreement between the signal profiles provided by the model and the MD simulation. Before fitting the model to the MD signals, we demonstrate a general aspect of 2D THz-Raman signals using the SL BO model and the optical Liouville paths in a simple case. This serves as a guid to subsequent analysis.

### iv.1 General aspects of 2D THz-Raman signals

Here, we elucidate several aspects of the signal components in 2D THz-Raman spectroscopy using the SL BO model and optical Liouville paths. As explained in Sec. I, we can express the dipole and polarizability in terms of the collective coordinate as and , respectively. If the inter-molecular modes are both Raman and THz active, the three THz-Raman signals are expressed as

(31) | ||||

(32) |

and

(33) |

with , where the anharmonic component is expressed as

(34) |

Similarly the nonlinear polarizability components are given by

(35) | ||||

(36) |

and

(37) |

In the harmonic LL BO case, the above THz-Raman signals can be calculated analytically, and we have , ,

(38) |

and

(39) |

where is the first-order response function of the harmonic BO system. For an isolated oscillator with frequency , we have . As explained in Sec. I, the term vanishes in the harmonic case because there is an odd number of Gaussian integrals involved in the response function.Okumura-JCP-1996-105 (); Okumura-CPL-1997-277 () Moreover, the term vanishes because of the cancelation of possible optical Liouville paths, as explained in Appendix B. While becomes large for large anharmonicity ,Okumura-JCP-1996-105 (); Okumura-CPL-1997-277 () the contribution of remains small due to this cancelation. Thus, in the anharmonic case, we can estimate from the RTT measurement. Then, by subtracting from the TRT and TTR signals, and , we can evaluate and separately. Because each contribution arises from the corresponding optical process, we can elucidate the key features of the inter-molecular motion from them. In the 2D Raman case, contrastingly, because all of the contributions appear together in a single observable, we cannot carry out such analysis.

To more clearly elucidate the characteristic of the 2D THz-Raman signals, in the figures, we display the 2D profiles of the , , and components separately in the slow modulation case ( and ) and the fast modulation case ( and ), as obtained from the SL BO model (, ) with and . Note that, because the effective system-bath coupling strength becomes weaker as the modulation becomes faster, we change both and to elucidate a pure non-Markovian effect.Tanimura-JCP-accepted () Although the contribution is minor, we plot the component in Appendix B. The relative intensities of each component evaluated in the harmonic and anharmonic SL BO cases are presented in Table. 1. We then analyze each profile using the optical Liouville paths (the double-sided Feynman diagrams).

Potential | Modulation | ||||
---|---|---|---|---|---|

harmonic | slow | — | |||

harmonic | fast | — | |||

anharmonic | slow | ||||

anharmonic | fast |

The 2D profile of the RTT component is presented in Appendix B.

#### iv.1.1 The TRT component,

In Fig. 2, we present the component for (a) the slow and (b) the fast modulation cases calculated using the harmonic SL BO model (). It is seen that the peak profiles are symmetric along the line, as can be deduced from the analytical expression for the LL case, Eq. (38). The peaks in the slow modulation case stretch in the direction, whereas those in the fast modulation case stretch in the direction, where is the fundamental frequency and an integer. In the slow modulation case, there are elongated peaks called ”echo peaks” along the direction.

Although our simulation results are fully classical, these profiles can be interpreted easily using the quantum Liouville paths for the TRT process depicted in Fig. 3. There, an energy eigenstate of the harmonic potential is denoted by , and we have depicted cases starting from the vibrational ground state as examples. The dipole operator , represented by the red circles in the diagram, converts the state into and through single quantum (SQ) excitations, while the nonlinear polarizability operator , represented by the blue double circles, converts the state into and through double quantum (DQ) excitations or maintains the same state through the zero quantum (ZQ) excitation. Because the final state, appearing after the last laser interaction, must be a population state , due to the trace operation involved in the response function, the possible processes are limited to the cases depicted in Figs. 3(i)-(iii) and their conjugate diagrams. While the double circle in the diagram (i) involves the transition , that in the diagram (ii) involves the transition . Thus, although these paths have same phases with opposite signs, they do not cancel. These six components constitute the signal expressed in Eq. (38) in the isolated oscillator case.

The phases for the paths (i) and (ii) in Fig. 3(i), expressed as , are the same because the net transition for these paths is a ZQ transition, while that for path (iii), expressed as , is different because the net transition for this path is a DQ transition. In the slow modulation case, the signal in Fig. 3 (iii) can rephase and become strong along the direction, but in the fast modulation case, coherence is lost due to rapid changes in the fundamental frequency, as illustrated in Fig. 1, and the signal decays quickly. Thus, we observe the echo signal in the TRT component in the slow modulation case, whereas we observe the chain of peaks along the direction in the fast modulation case, due to the processes depicted in Figs. 3 (i) and (ii).

#### iv.1.2 The TTR component,

We now discuss the component, presented in Fig. 4. The prominent features of this signal are the appearance of elongated peaks parallel to the axis and the appearance of peaks along the direction for small values of . The optical Liouville paths for are presented in Fig. 5. Three other paths can be obtained by exchanging the left and right arrows. These six components constitute the signal expressed in Eq. (39) in the isolated oscillator case.

We observe the population states during the period of Figs. 5(i) and 5 (ii) and the coherent states created by the two SQ excitations during the period of Figs. 5(iii). While all of the diagrams in Fig. 5 exhibit the oscillation during the period, only that in Fig. 5(iii) exhibits the oscillation during the period. Due to the SL interaction, the high-frequency oscillation of the coherent state appearing in Fig. 5(iii) decays quickly, while the population state appearing in Fig. 5(ii) remains for a long time. Thus, we observe the peaks along the direction for a short period, whereas elongated peaks appear along the direction.

#### iv.1.3 The AH component,

The contribution of the component arises only in the anharmonic case. In Fig. 6, we display the signal for in the (a) slow and (b) fast modulation cases, respectively. Although the definition of the response function is different, the optical Liouville paths for are similar to those for and . To explain the reason for this, we consider energy eigenstates in the case of an anharmonic potential with eigenenergy . The fundamental frequency between and is denoted by , whereas that between and is denoted by . In addition to such frequency shifts, the anharmonicity changes the roles of the dipole and polarizability operators. While the component involves only and , they can induce ZQ and DQ transitions, in addition to the SQ transition, because , , and are all nonzero in the anharmonic case. Because both THz and Raman laser pulses can induce ZQ and DQ transitions, the diagrams for include all of the diagrams in Figs. 3 and 5, while the resonant frequency for the period becomes , and that for the period is , , , or . The rephasing processes depicted in Fig. 3 (iii) occur only rarely, even in the slow modulation case, due to the anharmonicity, and this contribution can therefore be ignored. For this reason, the 2D profiles of the component are similar to those of and , but without the echo peaks described by Fig. 3 (iii). As in the TTR case, the contribution from the diagram depicted in Fig. 5 (iii) decays quickly. Thus, we observe elongated peaks in the direction that arise from population decay during the period.

### iv.2 The MD and fitted results

In order to concretely study hydrogen-bonding dynamics in 2D THz-Raman spectroscopy, we selected formamide, water, and methanol. While water and methanol exhibit high-frequency librational motion arising from hydrogen bonding, formamide exhibits only low-frequency inter-molecular motion. For each liquid, we explore the parameter values of the LL+SL BO model to fit the RTT, TRT and TTR signals to the MD simulation results. We fixed the anharmonicity to , while we varied the nonlinear polarizability . This was done because the intensity of the component is proportional to ,Okumura-JCP-1996-105 (); Okumura-CPL-1997-277 () and only the ratio of and is important to elucidate the roles of anharmonicity and non-linear polarizability. Because the coupling strengths of the LL and SL interactions are determined by and , respectively, we also fixed and varied and in the fit. The best fits of these parameters are presented in Table 2. Note that because the SL interaction increases the effective system frequency, Steffen-JPSJ-2000-69 (); Tanimura-JPSJ-2000-69 () the fitted do not correspond to the resonant frequency estimated from the 1D IR results obtained from the MD simulation.

Molecule | |||||
---|---|---|---|---|---|

Formamide | |||||

Water | |||||

Methanol |

#### iv.2.1 Formamide

In Fig. 7, we display (a) the MD results and (b) the fitted results from the LL+SL BO model for the (i) RTT, (ii) TRT, and (iii) TTR cases of liquid formamide. It is seen that the LL+SL BO results capture the essential features of the 2D THz-Raman spectra. The characteristic feature of the MD signals for formamide is the elongation of the peaks along the axis observed in Figs. 7(a-i) and (a-ii). By adapting the analysis applied to and , we conclude that this elongation arises from the slow population decay during the period described by the diagrams in Figs. 3 (i), 3 (ii), 5 (i), and 5 (ii). The echo peaks do not appear in Figs. 7(a-ii) and (b-ii), because the modulation in this case is so fast that dephasing cannot take place. In the TTR case depicted in Figs. 7(a-iii) and (b-iii), the peaks along the direction appear for small values of and due to the contribution of the diagram presented in Fig. 5 (iii). The similarity between the MD and model results indicates that the collective modes are subject to fast frequency modulation, as illustrated in Fig. 1(a).

#### iv.2.2 Water and Methanol

In Figs. 8 and 9, we display (a) the MD results and (b) the LL+SL BO results for the (i) RTT, (ii) TRT, and (iii) TTR cases of liquid water and methanol, respectively. We reproduced all of the 2D profiles obtained from the MD simulation for these liquids using the LL+SL BO model with the parameter values listed in Table 2. In these liquids, the vibrational modes near are the librational modes.Bosma-JCP-1993-98 (); Larsen-JCP-2006-125 ()

The characteristic features of the MD signals for these two liquids are the elongation of the peaks along the axis in the RTT and TRT signals and the echo peaks along the line in the TRT signals. The existence of these peaks indicates that the collective modes are inhomogeneously distributed, as depicted in Fig. 1 (b), due to slow frequency modulation. Because water and methanol exhibit high-frequency librational motion caused by hydrogen bonds, we observe an oscillatory feature in the 2D signal. This contrasts with the formamide signal, which decays quickly.

## V Conclusion

Using a Brownian oscillator (BO) model with a linear-linear (LL) and square-linear (SL) system-bath interaction, we analyzed the 2D RTT, TRT, and TTR signals for formamide, water, and methanol liquids obtained from full MD simulations. The classical hierarchal equations of motion approach was used to calculate the 2D signals with the LL+SL BO model under non-perturbative and non-Markovian conditions. By fitting the anharmonicity of the potential, the LL and SL coupling strength, and the inverse noise correlation time of the LL+SL BO model to the results of the MD simulations, we reproduced each of the simulated signal profiles by capturing their characteristic features. We found that the profile of formamide liquid can be accounted for by homogeneously distributed oscillators, whereas the profiles of water and methanol liquids can be accounted for by inhomogeneously distributed oscillators. Due to hydrogen bond interactions, water and methanol exhibit oscillating echo peaks. We were able to describe the two cases as the cases of fast and slow modulation of the anharmonic LL+SL BO model. The key feature of the present model that allows it to account for the simulated signal is the existence of the SL interaction. We were able to describe the complex inter-molecular motion with this simple model because it is capable of describing the collective modes under conditions ranging from homogeneous to inhomogeneous in a unified manner through variation of the noise correlation time. The success of the present study indicates that the LL+SL BO model captures the essence of the inter-molecular motion.

Finally, we briefly discuss some extensions of the present study. First, note that there does exist some discrepancy between the MD simulation results and the LL+SL BO model results. It should not be too difficult to decrease this discrepancy. For example, in the case of water, we should be able to improve the description of the signal profiles by introducing the second mode, which is , in the BO model. The significance of the fitted frequencies, , should be clarified through the normal mode analysis of MD simulations. Secondly, as shown in previous studies, the LL+SL BO model can be applied not only to 2D THz-Raman spectroscopy, but also to 2D IR spectroscopy for both inter-molecularTanimura-ACR-2009-42 () and intramolecular vibrations.Ishizaki-JCP-2006-125 (); Ishizaki-JPC-2007-111 (); Sakurai-JPCA-2011-115 () An extension to bridge between the inter- and intra-molecular modes using an extension of the present model should also be possible. We leave such extensions to future studies, in accordance with progress in experiments and simulations.

###### Acknowledgements.

Financial support from a Grant-in-Aid for Scientific Research (A26248005) from the Japan Society for the Promotion of Science is acknowledged. T. Ikeda is thankful to T. Hasegawa and A. Sakurai for their suggestion to solve the hierarchal Fokker-Planck equations in the Hermite representation.## Appendix A The Hierarchal Fokker-Planck Equations in the Hermite Representation

In the case of a strong system-bath coupling strength , Eqs. (18) and (19) converge very slowly as difference equations with a discrete mesh in the phase space. By expanding in terms of Hermite functions in the momentum direction, Risken-Book () the convergence can be improved. In this case, the HEOM become simultaneous equations for the coefficients of the expression. We expand the distribution function as follows:

(40) |

where is the th Hermite function,

(41) |

with is the th Hermite polynomial and .

The Liouvillian of the system (given in Eq. (20)) and the relaxation operators (given in Eqs. (21) and (22)) are expressed as

(42) | ||||

(43) |

and

(44) |

where ,

(45) |

and

(46) |

Then, the equations of motion for the coefficients are reduced to

(47) |

for and

(48) |

To carry out the numerical calculation, we chose so as to satisfy and solved Eqs. (47)-(48) as simultaneous equations.