Two-dimensional (2D) vibrational spectroscopy has emerged as one of the most important experimental techniques useful to study the molecular structure and dynamics in condensed phases. Theory and computation have also played essential and integral roles in its development through the nonlinear optical response theory and computational methods such as molecular dynamics (MD) simulations and electronic structure calculations. In this article, we present the fundamental theory of coherent 2D vibrational spectroscopy and describe computational approaches to simulate the 2D vibrational spectra. The classical approximation to the quantum mechanical nonlinear response function is invoked from the outset. It is shown that the third-order response function can be evaluated in that classical limit by using equilibrium or non-equilibrium MD simulation trajectories. Another simulation method is based on the assumptions that the molecular vibrations can still be described quantum mechanically and that the relevant molecular response functions are evaluated by the numerical integration of the Schrödinger equation. A few application examples are presented to help the researchers in this and related areas to understand the fundamental principles and to use these methods for their studies with 2D vibrational spectroscopic techniques. In summary, this exposition provides an overview of current theoretical efforts to understand the 2D vibrational spectra and an outlook for future developments.
I. INTRODUCTION
Resonance is a phenomenon ubiquitous in nature. Molecular spectroscopy is an experimental tool that utilizes one of the resonance phenomena that involves an interaction of oscillating charged particles in a given molecule with an external electromagnetic field whose frequency (ωwave) is close to that (ωmolecule) of a vibrational or electronic oscillation determined by the associated interaction potential between constituent charged particles in a given polyatomic molecule in condensed phases. Such resonant field-matter interaction is manifested in the complex quantum transition amplitude that is approximately proportional to (ωwave − ωmolecule + iγmolecule)−1, where γmolecule determines the linewidth of such a transition band that reflects complicated molecular interactions with fluctuating internal and bath degrees of freedom. As the wave frequency ωwave approaches the molecular oscillation frequency ωmolecule, the corresponding quantum transition amplitude increases, which has been referred to as a quantum resonant field-matter interaction process.
Vibrational spectroscopy, one of the molecular spectroscopic techniques, involves quantum transitions of vibrational degrees of freedom. The most widely used infrared (IR) spectroscopy measures the vibrational transition amplitude of normal modes via their resonant interaction with infrared (e.g., near IR, mid-IR, far-IR, and THz) radiation. The vibrational Raman spectroscopy is a resonance-enhanced scattering of the incident electronically non-resonant field when the beat frequency between the incident radiation field and the inelastic Raman scattering field is close to molecular vibrational frequency. Although IR-vis sum-frequency-generation (SFG) spectroscopy for molecular systems on the surface or at the interface with broken centrosymmetry is a nonlinear optical (three-wave-mixing) spectroscopy, it still involves a single vibrational transition whose amplitude is essentially determined by the quantum correlation between the transition dipole and the transition polarizability of a given vibrational mode.
Coherent multidimensional vibrational spectroscopy1–5 is an extension of these linear spectroscopic methods and involves multiple vibrational transitions that are temporally separated. To observe and analyze the correlation between the distinctive vibrational transitions, one usually applies a series of coherent laser pulses with specific phase relations to induce such transitions. The multiple transitions are detected and presented as a function of corresponding frequency variables, making the spectrum multidimensional. In particular, coherent two-dimensional (2D) vibrational spectroscopy employs three femtosecond laser pulses in the infrared (IR) or visible frequency range to induce the third-order polarization in molecular systems. The signal due to the polarization is measured and presented in two frequency dimensions conjugates to the time intervals between the first and second pulses (τ) and between the third pulse and the detection (t), respectively. A series of 2D spectra at different time intervals between the second and the third pulse (T) reveals the molecular dynamics (MD) taking place during that time interval T that is called the waiting or population time.
2D vibrational spectroscopy provides ultrafast time resolution along the waiting time (T) axis as well as high spectral resolution of individual peaks along the diagonal in the 2D frequency space and their couplings in the form of off-diagonal cross peaks. Therefore, the method conveys rich information on molecular systems such as homogeneous (anti-diagonal) and inhomogeneous (diagonal) spectral broadening effects, vibrational anharmonicity, spectral diffusion,6 and intermode coupling strength and its temporal variation. Over the past two decades, 2D vibrational spectroscopy has been extensively used to study the structure and dynamics of small peptides, proteins, DNA, and lipid bilayers, energy transfer dynamics, hydrogen-bonding (H-bonding) structure and dynamics of liquid water, and its isotopologues, and configurational and H-bonding dynamics of biomolecules.
For coupled multi-oscillator systems, vibrational couplings through space via intermolecular interactions or through bond via anharmonicities in the multi-dimensional potential energy surface are crucial in understanding vibrational dynamics.9,10 However, since the linear vibrational spectrum is mainly determined by the harmonic properties, e.g., vibrational frequency and transition dipole moment, of the oscillators, it is difficult to quantitatively extract such weak features, e.g., vibrational coupling constants and potential anharmonic coefficients of coupled oscillators in condensed phases, from those linear spectra.8 On the other hand, the coherent multi-dimensional vibrational spectroscopic methods have been found to be of exceptional use because they enable one to explore these mode-mode couplings and provide crucial information on the structural dynamics through the time-dependent changes of the spatial proximity and relative orientation of chromophores. Note that the nonlinear vibrational response functions determining amplitudes and lineshapes of coherent multidimensional vibrational spectra vanish for perfect harmonic oscillators.9,11–14 One of the most popular techniques is the 2D IR spectroscopy,1 which is a four-wave-mixing method capable of providing information on the molecular nonlinear response function that is given by multi-time correlation functions of transition dipole moments. Again, using the linear expansion form of electric dipole moment with respect to vibrational coordinates and employing perturbation theory treating electric and mechanical anharmonicities as weak perturbations to harmonic oscillators, it was possible to rewrite the nonlinear response functions in terms of vibrational Green functions and such perturbation terms.11,12,15,16 This method was known as a linked diagram theory for the multidimensional vibrational response function. Although such theoretical approaches were popular in the early time of theory on multidimensional vibrational spectroscopy, the intrinsic difficulties in describing the anharmonicity-induced frequency shift of the excited state absorption contribution to the 2D IR spectrum, for example, and in accurately calculating vibrational (both mechanical and electrical) anharmonicities for oscillators in solutions, prohibited its wide use later.
Although 2D IR spectroscopy employing three incident IR laser pulses is one of the most widely used forms of coherent multidimensional vibrational spectroscopy, different varieties of the method have also been developed and applied for specific purposes. Examples include surface specific 2D SFG spectroscopy,17,18 2D Raman19 and terahertz spectroscopy,20–24 2D IR-IR-visible spectroscopy,3,9,13,25–29 and so on. In this article, because the basic theories for these variants are similar to one another, we mainly focus on the theory and computation of 2D IR spectroscopy and reference is made to these variants when appropriate.
By necessity, the exposition in this paper cannot be exhaustive. However, it is hoped to serve as an up-to-date introduction and a basis for elaborate research especially with the aid of references provided. In Sec. II, we briefly introduce the third-order response function formalism, which has the central importance in the description of 2D vibrational spectroscopy. In Sec. III, classical mechanical approaches to calculate 2D vibrational spectra are presented based on the classical limit of the quantum vibrational third-order response function. In Sec. IV, we present a quantum mechanical (QM) method to simulate the spectra by numerical integration of the Schrödinger equation (NISE). Finally, we conclude with perspective on future development.
II. FUNDAMENTALS
In general, spectroscopic measurements are conducted in two steps: the preparation step where the molecular system is excited by incident radiations and the detection step where the generated signal is measured and presented. In 2D vibrational spectroscopy, the system is irradiated with three coherent laser pulses, usually in the infrared (IR) frequency region, in the preparation step and then the generated signal field is detected and presented in two frequency dimensions, representing two distinct coherence oscillations separated by a waiting (population) time. It is a kind of four-wave-mixing spectroscopy due to the fact that the signal field arises from three preceding field-matter interactions that are each linear in the field. In each of the four field-matter interaction events, a quantum transition takes place between vibrational states of the system, on either the ket or the bra side of the system density operator (see below). Depending on the configuration of the optical laser pulses such as the frequency, the direction of propagation (wave vector), and polarization, as well as the detection methods, different quantum transition pathways can be selectively generated and measured.8
The spectroscopic observables are determined by the third-order optical response function of the molecular system that includes all distinctive quantum transition sequences consistent with the incident radiations and the detection method. The response function naturally emerges from the quantum mechanical time-dependent perturbation theory treatment of the molecular system in the presence of the three perturbative light-matter interactions of the preparation step. In this section, we sketch the theoretical procedure and present key results relevant to 2D vibrational spectroscopy.
A. Third-order optical response function
B. Response function of a three-level system
2D vibrational spectroscopy usually induces transitions up to the second vibrational excited state. Therefore, a three-level system with eigenstates , , and is a useful model for the response function relevant to 2D vibrational spectroscopy. Because the third-order response function vanishes for a harmonic oscillator, the model system must represent an anharmonic oscillator where the fundamental transition frequency ωeg is slightly larger than ωfe.
Double-sided Feynman diagrams illustrating the six contributions of Eq. (12). Time is running from bottom to the top with the vertical lines illustrating the ket and the bra wave functions. Interactions with applied laser fields are illustrated with full arrows, while emitted light is illustrated with dashed arrows. The numbers (zero, one, and two) indicate the number of vibrations excited. The response functions R2A and R4 are identified as ground state bleach as the system is in the ground state during the waiting time t2, while R3 and R1A are stimulated emission contributions, and R1B and R2B are excited state absorption contributions.
Double-sided Feynman diagrams illustrating the six contributions of Eq. (12). Time is running from bottom to the top with the vertical lines illustrating the ket and the bra wave functions. Interactions with applied laser fields are illustrated with full arrows, while emitted light is illustrated with dashed arrows. The numbers (zero, one, and two) indicate the number of vibrations excited. The response functions R2A and R4 are identified as ground state bleach as the system is in the ground state during the waiting time t2, while R3 and R1A are stimulated emission contributions, and R1B and R2B are excited state absorption contributions.
This general formulation for three-level systems is very powerful and can be used to understand the effect of dynamics on the two-dimensional lineshapes,36 which, for example, allows the distinction between homogeneous and inhomogeneous line-broadenings and defines rules for extracting the correlation functions describing the bath dynamics in the limit of the Gaussian dynamics, for example, using the nodal or center line slope.6,37 This is illustrated in Fig. 2.
Illustration of the 2D IR spectra with slow (left) and fast (right) dynamics. The corresponding linear spectra, which have little sensitivity to dynamics, are shown on the top. The diagonal peaks shown with red contours are bleaching signals (R1A, R2A, R3, and R4), while the peaks with blue contours, which are shifted by the anharmonicity below the diagonal, arise from the absorption of the excited state (R1B and R2B). When the bath dynamics scramble the vibrational frequencies faster than the waiting time, all correlation between excitation and detection frequencies is lost and round peaks are obtained, while elongated peaks arise from residual memory.
Illustration of the 2D IR spectra with slow (left) and fast (right) dynamics. The corresponding linear spectra, which have little sensitivity to dynamics, are shown on the top. The diagonal peaks shown with red contours are bleaching signals (R1A, R2A, R3, and R4), while the peaks with blue contours, which are shifted by the anharmonicity below the diagonal, arise from the absorption of the excited state (R1B and R2B). When the bath dynamics scramble the vibrational frequencies faster than the waiting time, all correlation between excitation and detection frequencies is lost and round peaks are obtained, while elongated peaks arise from residual memory.
III. CLASSICAL MECHANICAL APPROACHES TO NONLINEAR RESPONSE FUNCTIONS
In Sec. II, we derived the quantum mechanical linear and nonlinear response functions. Although these equations are formally exact, fully quantum mechanical simulations are still impractical for systems with many degrees of freedom even with state-of-the-art supercomputers. Thus, applicable and efficient methods are required for the calculation of nonlinear response functions for such systems. In this section, we first summarize the classical mechanical approaches to the calculation of the linear and nonlinear response functions, which can be expressed as (multi-)time correlation functions on equilibrium and nonequilibrium trajectories, and then present some results of the 2D IR and pump-probe spectra obtained from the approaches.
A. Equilibrium molecular dynamics approach
Jeon and Cho have investigated the 2D IR spectra of intramolecular vibrations in water.49,54 In their studies, the quantum mechanical/molecular mechanical (QM/MM) approach has been employed for the accurate description of the intramolecular vibrations of solute molecules: deuterated N-methylacetamide (d-NMA) and HOD molecules have been described with the semiempirical PM3 method and scc-DFTB potential, respectively. The flexible TIP3P and SPC/Fw models have been used for solvent water molecules, respectively. Figure 3 shows the 2D IR spectra of the OD stretch of the HOD molecule in a water cluster. In these calculations,54 Jeon and Cho have improved the computational efficiency by exploiting the time reversibility of trajectory and have successfully calculated the 2D IR spectra. The positive and negative peaks are found at (∼2650 cm−1, ∼2700 cm−1) and (∼2600 cm−1, ∼2550 cm−1), respectively. The positive peak corresponds to the stimulated emission and the ground state bleaching of the amide I band, whereas the negative peak is related to the excited state absorption. Note that various general features of the experimental 2D IR spectra are reproduced in the calculated 2D spectra based on classical mechanics: the vertical splitting of positive and negative peaks due to the vibrational anharmonicity, the diagonal elongation of the signal reflecting the inhomogeneity of the solvation environment, the change in lineshape with waiting time due to the spectral diffusion, and the decrease in the tilt angle of the nodal line. The classical 2D spectra also show that the main decay component of the nodal line with a ∼1.6 ps time scale obtained from the calculated 2D spectra is comparable to experimental results.55
Absorptive 2D IR spectra of the OD stretch of the hydrated HOD molecule at different waiting times T. All figures are drawn on the same amplitude scale. Diminishing tilt angles of the nodal line and the signal amplitudes can be noticed with increasing waiting time T, indicating the loss of frequency memory and the vibrational population relaxation, respectively. Reprinted with permission from J. Jeon and M. Cho, J. Phys. Chem. B 118, 8148 (2014). Copyright 2014 American Chemical Society.
Absorptive 2D IR spectra of the OD stretch of the hydrated HOD molecule at different waiting times T. All figures are drawn on the same amplitude scale. Diminishing tilt angles of the nodal line and the signal amplitudes can be noticed with increasing waiting time T, indicating the loss of frequency memory and the vibrational population relaxation, respectively. Reprinted with permission from J. Jeon and M. Cho, J. Phys. Chem. B 118, 8148 (2014). Copyright 2014 American Chemical Society.
B. Nonequilibrium molecular dynamics approach
1. Nonequilibrium finite-field method
In Subsection III A, we described the prescription for the linear and nonlinear response functions based on the equilibrium molecular dynamics approach. Since the stability matrix is required in the approach, high computational costs and numerical instability would make it difficult to calculate nonlinear response functions, when many degrees of freedom have to be considered explicitly, for example, the low frequency intermolecular motions and the OH stretches in liquid H2O. Thus, other approaches circumventing the use of the stability matrix would be required for such systems. In this subsection, we introduce the nonequilibrium molecular dynamics approach in which the response functions are evaluated by considering (sequential) external fields, as in real experiments.
2. Hybrid approach combining equilibrium and nonequilibrium finite-field methods
The nonequilibrium finite-field method does not need the stability matrix for the calculation of nonlinear responses. However, the approach still demands high computational costs; i.e., extra three and seven nonequilibrium trajectories as well as the equilibrium trajectory are required for the second- and third-order nonlinear response functions, respectively.
It should be noted that the experimental third-order nonlinear spectra are obtained under different phase matching conditions [see Eq. (8)]. For example, the 2D IR spectra consist of the responses emitted along the wave vector directions kI = −k1 + k2 + k3 and kII = +k1 − k2 + k3. However, the calculated third-order response function consists of all the wave vector components generated by three incident electric fields, i.e., kIII = +k1 + k2 − k3 and kIV = +k1 + k2 + k3 as well as kI and kII. Therefore, the components with kIII and kIV have to be eliminated from the calculated third-order responses for the 2D IR spectrum. The responses arising from the individual phase matching conditions can easily be obtained by using the nonequilibrium molecular dynamics approach in which external electric fields are explicitly considered. Thus, various third-order nonlinear spectroscopies, such as the three-pulse photon echo peak shift and the pump-probe spectrum, can also be calculated.61 In the earlier studies, kIII and kIV components have been eliminated by exploiting the fact that the three-dimensional Fourier transformed spectra of kIII and kIV components show higher frequency oscillation in ω2 than those of kI and kII.47–49,54
Hasegawa and Tanimura have applied the backward-forward sampling method to the 2D IR spectra of liquid HF.47 As shown in Fig. 4, the positive and negative peaks of intermolecular libration motion are located at (550 cm−1, 550 cm−1) and (550 cm−1, 400 cm−1), respectively. The frequency correlation of the libration motion of liquid HF was found to be lost with a time scale of ∼200 fs. It was also found that the width of the anti-diagonal line of the peak is narrow compared with that of liquid water, indicating the suppressed homogeneous broadening in liquid HF.
Correlation spectra of liquid HF at t2 = (a) 0 fs, (b) 20 fs, (c) 100 fs, (d) 200 fs, (e) 300 fs, and (f) 400 fs. The diagonal and off-diagonal peaks in each spectrum arise from fundamental and anharmonic oscillations, respectively, of the hydrogen bond network in liquid HF. Reprinted with permission from T. Hasegawa and Y. Tanimura, J. Chem. Phys. 128, 064511 (2008). Copyright 2008 AIP Publishing LLC.
Correlation spectra of liquid HF at t2 = (a) 0 fs, (b) 20 fs, (c) 100 fs, (d) 200 fs, (e) 300 fs, and (f) 400 fs. The diagonal and off-diagonal peaks in each spectrum arise from fundamental and anharmonic oscillations, respectively, of the hydrogen bond network in liquid HF. Reprinted with permission from T. Hasegawa and Y. Tanimura, J. Chem. Phys. 128, 064511 (2008). Copyright 2008 AIP Publishing LLC.
Yagasaki and Saito investigated the fluctuation and relaxation of the intermolecular motions in liquid water by using the hybrid method with the SPC/E model for water molecules.48 Figure 5 shows the 2D IR spectra of the intermolecular motions of liquid water at several waiting times. The intermolecular translational and librational motions are found below and above 300 cm−1, respectively. As shown in Fig. 5, the tilt angle of the nodal line between the positive and negative peaks of the librational motion decays with a time constant of ∼110 fs. Yagasaki and Saito found that the frequency fluctuation of librational motion becomes three times slower when the translational motions of individual molecules are frozen, indicating that the ultrafast frequency fluctuation of the librational motion arises from the coupling between the translational and librational motions. It was found that the three-pulse stimulated photon-echo signal of the librational motion in water decays with time scales of ∼20 and ∼100 fs, in which the slower time scale is similar to that obtained from the change of the nodal line tilt angle.48
2D IR correlation spectra of intermolecular motion of liquid water at several waiting times. The tilt angle is defined as the angle between the ν1 axis and the nodal line of the libration peak as shown in (a) and is found to decay with a time scale of 115 fs. Reprinted with permission from T. Yagasaki and S. Saito, J. Chem. Phys. 128, 154521 (2008). Copyright 2008 AIP Publishing LLC.
2D IR correlation spectra of intermolecular motion of liquid water at several waiting times. The tilt angle is defined as the angle between the ν1 axis and the nodal line of the libration peak as shown in (a) and is found to decay with a time scale of 115 fs. Reprinted with permission from T. Yagasaki and S. Saito, J. Chem. Phys. 128, 154521 (2008). Copyright 2008 AIP Publishing LLC.
The energy relaxation process of the librational modes has been examined by the pump-probe signals calculated from the third-order response functions.61 It was found that the absorption changes at 700 and 800 cm−1 are induced by the initial decrease of absorption due to the ground state bleach and stimulated emission followed by the fast energy relaxation of the librational motion to low frequency modes within ∼60 fs and the subsequent slow thermalization that corresponds to its energy relaxation to hot ground states occurs in 500 fs. These time scales of the two relaxation processes are consistent with the experimental results.62–64
As mentioned above, once third-order response functions with individual phase-matching conditions are obtained, any third-order nonlinear spectrum can be calculated in principle. However, despite these developments of efficient computational approaches, the calculation of the third-order response function is still expensive and computationally demanding. In this regard, it should be noted that other efficient computational methods to estimate mode-to-mode vibrational energy relaxation rates have been developed by Yagasaki and Saito65,66 and Jeon and Cho.67,68
Recently, the nonequilibrium molecular dynamics method has also been applied to the intramolecular vibrational modes of H2O in water.69 In this study, the 2D IR spectra of the OH stretch and the HOH bending modes were calculated with the TTM3-F interaction potential, which is the ab initio-based charge transferable, flexible, and polarizable Thole-type model to describe the inter- and intra-molecular interactions in liquid water. The waiting time-dependent changes of the OH stretch spectra indicate the heterogeneous dynamics of the local H-bonding network which has also been found in an experimental result.70 The spectrum loses diagonal correlation at waiting times larger than 200 fs, which is consistent with the other quantum mechanical calculation results.71,72 The 2D IR spectra of the HOH bend are found to have a smaller nodal line tilt angle at T = 0 fs compared to the OH stretch spectra and become parallel to the ω1 axis within ∼400 fs.
The pump-probe spectra of the OH stretch and the HOH bend in water have also been calculated with the TTM3-F interaction potential.69 The positive and negative peaks of the pump-probe spectra of the OH stretch at (ωpump, ωprobe) = (3600 cm−1, 3225 cm−1) and (3600 cm−1, 3600 cm−1), corresponding to the 1 → 2 and 0 → 1 transitions, show the initial decay with a time constant of 240 fs followed by the slow relaxation to the hot-ground state. The calculated time scale of the initial decay is in agreement with experimental results of 200-270 fs.73–75 The pump-probe spectra of the HOH bend at (1650 cm−1, 1475 cm−1) and (1650 cm−1, 1650 cm−1) decay with a time constant of 250 fs, which is in good agreement with the experimental results of 170–260 fs.62,64 Furthermore, the frequency-resolved kinetic energy analysis revealed the detailed relaxation processes after the excitations of the OH stretch and HOH bend of H2O in water.
C. Limitations
In this section, the classical simulation methods for the nonlinear response functions were presented with their applications to the calculations of the 2D IR and pump-probe spectra of liquids. Before closing this section, we summarize the limitations of the classical treatment of nonlinear spectral simulations.
The validity of the classical 2D IR spectra has been investigated by several groups.76–81 It should be noted that the classical nonlinear response functions are not stable, for example, for integrable systems and systems without dissipation.77,78 Sakurai and Tanimura examined the quantum effects on the simulated IR and 2D IR spectra of a Morse oscillator in a harmonic bath by solving the quantum and classical hierarchical equations of motion (HEOM).80 They found that the classical 2D IR spectra are a good approximation of the quantum 2D IR spectra when the system (vibration) is largely modulated by the bath, i.e., via a strong system-bath coupling or a fast bath modulation even in a weak system-bath coupling. On the other hand, there are significant differences between the quantum and classical mechanically calculated spectra when the modulation due to the system-bath interactions is weak or slow. Very recently, Reppert and Brumer have also shown that the classical 2D IR spectra of a Morse oscillator, which mimics the amide I vibration, can reproduce the qualitative features of the quantum 2D IR spectra very well.81 As shown in these studies, the validity of the classical mechanically calculated spectra depends on the system-bath coupling. Notably, in their study, it was shown that the anharmonicity of the classical signal is determined almost entirely by the frequency resolution afforded by the simulated scan time t3.
Classical nonlinear spectral simulations can reproduce spectra calculated with quantum methods qualitatively well when the system vibration is largely modulated by the bath, such as those examples discussed in this section. It should be noted, however, that any vibrational spectra are described as transitions from one vibrational level to another in quantum mechanics, whereas those are obtained from the fluctuation of dipole moment in classical mechanics. In quantum mechanics, wave functions or density matrices are determined by non-local information; i.e., the kinetic energy operator is expressed as the second derivative of coordinates in the coordinate representation. The anharmonic shift in quantum 2D IR spectra is thus determined by information about the whole potential energy profile. On the other hand, no discrete vibrational levels are considered in any classical approaches. In classical mechanics, a trajectory is determined by local information on coordinates and momenta. Consequently, the anharmonic shift corresponding to the frequency difference between the positive and negative peaks in the classical 2D IR spectra arises from the difference between the curvatures of trajectories perturbed by one and two electric fields. As a result, in many examples including those shown in this section, the anharmonic shifts in the classical 2D IR spectra tend to be smaller than those in the experimental or quantum mechanically calculated 2D IR spectra. In addition to the potential profile, the transition dipole moment is also important quantity for quantitative simulations of the IR and 2D IR spectra. In quantum mechanics, a spectral intensity is related to the transition dipole moments between discrete vibrational levels. On the other hand, the classical spectra are obtained by the dipole moment induced by vibrational and conformational changes. Accurate model potential and dipole moment are therefore extremely important to accurately simulate the 2D IR spectra of complicated molecular systems.
A fully quantum mechanical simulation of nuclear degrees of freedom including all the effects of surrounding thermal bath remains unpromising, despite the technological advancement of computers. Therefore, semiclassical approaches, such as (linearized) semiclassical-initial value representation,82 centroid molecular dynamics,83 and ring-polymer molecular dynamics,84 have been developed to accurately calculate the static and dynamic properties of molecules in condensed phases; for example, the IR spectrum of liquid water.85 A formalism for two-time correlation functions has also been developed.86 It is hoped that these semiclassical approaches will be able to provide a novel framework to calculate the 2D IR spectra incorporating the nuclear quantum effects by developing the formalism of nonlinear three-time correlation functions.87
IV. NUMERICAL INTEGRATION OF THE SCHRÖDINGER EQUATION
In practical calculations, the Hamiltonian in Eq. (30) only includes a change in the number of vibrations in the terms involving the coupling with the external electric field(s) and thermalization is not included. Therefore, the remaining part of the Hamiltonian is block diagonal. This allows treating the ground state, single excited states, and double excited states, separately. As there is only one ground state and the energy of this is set to zero, the time-evolution operator is the unit operator. When N vibrational degrees of freedom are treated, there are N single excited states and N(N + 1)/2 double excited states. The time-evolution in each excitation manifold can be evaluated independently in the corresponding harmonic basis, which, in principle, requires the evaluation of matrix exponentials of the corresponding dimensionalities, which is typically done by diagonalizing the Hamiltonian, multiplying the eigenvalues with −iΔt/ℏ, taking the exponent and transforming back to the original basis. However, in practice, the time-evolution matrix for the double excited states can be made more efficient as the time required for the direct evaluation of the time-evolution operator scales as N,6 which is a costly procedure to repeat even for moderately sized systems. Among the solutions to this problem are the nonlinear exciton propagation scheme,91 which utilizes an anharmonic perturbation of the harmonic solution, and an approach based on the Trotter algorithm making use of the sparse nature of the Hamiltonian of the doubly excited states.92,93 The latter allowed the calculation of the 2D IR spectra of systems up to 864 coupled vibrations for water ices94,95 as well as applications to the amide I band of full proteins.96–98 Numerous applications of the NISE method for 2D vibrational spectroscopy have been implemented when only one vibration is involved.99–102 A few different implementations applicable to coupled systems have been presented.88–93,103
A. Hamiltonian parameterization, frequency maps
The terms in the Hamiltonian in Eq. (30) have different origins, and their bath dependence can be parameterized in different ways. In the following, the different terms will be considered one by one. First, the vibrational frequency of a given mode depends on the local environment. The well-known Stark shift is a simple way to describe such frequency shift and fluctuation. Local electric fields modulate the vibrational frequency, and this type of effect can be parameterized by ab initio calculations of the vibrational frequency in different point charge environments of inhomogeneous electric fields. In this way, the solvent dependence of the vibrational frequencies of a number of common vibrations has been parameterized expressing the vibrational frequency in an expansion of the electric potential104–108 or electric field102,109–125 and sometimes the electric field gradient126–130 on specific points inside a molecule. The parameterization of the solvatochromic vibrational frequency shift parameters can either be obtained from ab initio calculations or empirical fitting to experimental data. These mappings then allow obtaining a frequency trajectory from a molecular dynamics trajectory. As typical frequencies fluctuate on a sub-picosecond time scale, the trajectories need to be sampled and saved every ∼20 fs and typical lengths are on the order of nanoseconds.
In reality, the solvatochromic vibrational frequency shift is not just an electrostatic effect. Charge transfer, dispersion forces, Pauli repulsion, polarization, and multipole effects may further play a role.121,125,131–137 Empirical models may compensate for this as long as other effects are correlated with electrostatic potential/field fluctuations. Ab initio calculations based on clusters, where the electrons of the solvent molecules are treated explicitly, may also be able to account for such effects in an efficient way. Detailed analysis of the amide I vibrations has shown that the various other effects tend to cancel out.132,136 For other systems, this may not be the case.
The anharmonicity also depends on the solvent, and this dependence has been included in mappings.102,126,128 Still, for many systems such as the amide I vibrations in proteins, the variation in the anharmonicity is small relative to the variation in the frequencies and the solvent dependence can be neglected using a constant anharmonicity value for each type of vibrational mode.
The couplings between different local modes, for long-range interactions, are well approximated by the transition-dipole coupling model.138–143 Transition charge coupling models have been developed for intermediate range couplings.130,144,145 These models account for some multipole contributions to the electrostatic interaction. For vibrational modes directly connected through covalent bonds, ab initio based coupling models have been developed. This has both been done for the OH stretch vibrations in the water molecule,146 for neighboring amide I modes in a peptide chain,117,118,130,145,147,148 and amide I and II modes in the same and neighboring peptide units of proteins.129,149 Ab initio modeling of the vibrational couplings in small peptides and DNA has also been reported based on the polarizable continuum model (PCM) description of solvent.150,151 For the peptide units, this coupling depends on the local configuration, as characterized by the Ramachandran angles. The developed mappings were thus made by varying these angles systematically for small di- or tri-peptides and then calculate the vibrational coupling through ab initio calculations. Just as for the vibrational frequencies, the couplings can be extracted from molecular dynamics trajectories using these mappings.
Figure 6 illustrates a comparison between the spectra simulated for the OH stretch of water molecules diluted in acetonitrile and experimental data. The evolution of cross-peaks between the high-frequency asymmetric OH-stretch and the low-frequency symmetric OH-stretch is well described by the theory. Furthermore, the difference in the dynamics of the two peaks, which was the focus of the study,152 is revealed by the center line slope analysis. From the simulations, this could be interpreted in terms of non-Gaussian dynamics as the bath dynamics for strong hydrogen bond environments is faster than for weak hydrogen bond environments. As strong hydrogen bonded OH-stretches absorb at lower frequencies, the fast dynamics contributes more to the low frequency symmetric vibration peak. While the anharmonicity of OH-stretch modes is around 200 cm−1, a sequential transition is observed with a much smaller apparent anharmonicity. This phenomenon arises as the large anharmonicity largely localizes the double excited states.
2D IR spectra for parallel laser polarization of the OH-stretch of water dissolved in acetonitrile at different waiting times t23 reproduced with permission from Jansen et al. J. Phys. Chem. A 113, 6260 (2009), Copyright 2009 American Chemical Society. (a) and (b) present experimental and simulated data, respectively. Dotted lines show the frequencies of symmetric and asymmetric stretching modes of the H2O molecule. The 2D spectra are normalized to the maximal amplitude. Red contours illustrate bleach signals, while the cyan contours illustrate excited state absorption. The contour lines are drawn at 10% steps of the maximal amplitude in each individual plot. Solid black illustrated the max lines.
2D IR spectra for parallel laser polarization of the OH-stretch of water dissolved in acetonitrile at different waiting times t23 reproduced with permission from Jansen et al. J. Phys. Chem. A 113, 6260 (2009), Copyright 2009 American Chemical Society. (a) and (b) present experimental and simulated data, respectively. Dotted lines show the frequencies of symmetric and asymmetric stretching modes of the H2O molecule. The 2D spectra are normalized to the maximal amplitude. Red contours illustrate bleach signals, while the cyan contours illustrate excited state absorption. The contour lines are drawn at 10% steps of the maximal amplitude in each individual plot. Solid black illustrated the max lines.
B. Mappings for interactions and the related techniques
A linear absorption spectrum can be calculated from the one-time (two-point) transition-dipole response function, while the 2D IR spectra are governed by a combination of six three-time (four-point) transition-dipole response functions related to rephasing and non-rephasing versions of the ground state bleach, stimulated emission, and excited state absorption pathways (Sec. II). The transition-dipoles of different states are reflected in the relative intensity of the peaks. Ultrafast reorientation of the transition-dipoles resulting from molecular reorientation153–156 and vibrational energy transfer72,153,157 can be monitored by measuring 2D IR spectra using different polarizations158,159 of the applied laser fields. Cross peak intensities in the 2D IR spectra depend on the relative orientation of the transition dipoles of the two involved resonances. This allows for the determination of molecular geometries160–162 and is powerful in the peak assignment. For many vibrations, the transition-dipole magnitude is fairly independent of the solvent environment. A clear exception of this is the OH-stretch vibration, which in water shows a five-fold change across the spectrum,163 with the very weak dipole strength of the free OH-stretch vibration and very strong absorption of a strongly H-bonded OH species. Therefore, it is crucial to include the variation of the transition-dipole in the modeling. This has been done very similar to how the solvent effects are accounted for through mapping for the vibrational frequencies.102,128,130,163
Raman spectra,102,146,164,165 sum-frequency generation,166–169 and 2D sum-frequency generation17,18,168,170–173 spectroscopy signals can be calculated using similar response functions as for linear absorption and 2D IR spectroscopies. The only major change is that one needs to replace transition-dipoles in the expressions with transition-polarizabilities. The possible polarization setups available increase as the transition-polarizabilities represent two interactions with external visible laser fields, while the transition-dipoles represent one interaction with an infrared laser field. This opens a possibility of directly probing molecular orientational distribution. Most importantly in the sum-frequency generation type experiments, which are only sensitive to non-centrosymmetric signal contributions,167,168,174,175 the orientation of the vibrations relative to the sample surface can be determined. Just as the transition-dipole, the transition-polarizability depends on the environment. This has been studied for the OH-stretch of water using an electric field mapping. Still, the solvent effect both for the OH-stretch and other systems is relatively small.102,146 The difference between the transition-dipole and the transition-polarizability solvent effects results in a non-coincidence effect147,176,177 making the infrared and Raman type spectra exhibit significantly different lineshapes. In particular, Raman and SFG spectra in the OH-stretch region are much more sensitive to free OH-species than infrared spectra.146
C. Limitations—Bandwidth, high-frequency modes
The quantum-classical methods discussed here have a number of important advantages and limitations. The methods are powerful in providing a quantum description of a given system and allow for treating essentially arbitrary bath dynamics, as provided by the combination of mappings and molecular dynamics simulations. For example, the methods allow describing non-Gaussian dynamics152,178–182 and chemical exchange101,155,179,183–192 arising when the bath coordinates are not harmonic. Commonly used methods invoking the second-order cumulant approximation34,193 or methods assuming a harmonic bath spectral density194–197 cannot account for such effects. The methods can use input from molecular dynamics simulations and are thus able to predict spectra starting from first principles. Sometimes it is not desirable to perform a full molecular dynamics simulation, and then stochastic models can be employed instead.183,198
The Hamiltonian in Eq. (30) does not allow for relaxations between the different excitation manifolds. This means that vibrational energy relaxation is intrinsically neglected. Furthermore, the methods account for the effect that the bath exerts on the system, while the feedback of the system to the bath when the system is on an excited state is not included.140,199–201 This results in a lack of correct thermalization in the quantum system, which leads to an equal population of all quantum states after equilibration.201 In practice, this means that artifacts can be expected at low temperatures, and when broad spectral ranges are considered, where the thermal energy is comparable to the width of the considered spectral region. In systems like water and ice, the OH-stretch band is broader than the thermal energy at room temperature, but still fairly good approximations of the spectra are found for such systems as long as the waiting time is not too long.72 For longer waiting times, persistent ground state bleach signals are observed in experiments resulting from the thermalization effect. A number of attempts to overcome this issue have been reported;140–143 however, including feedback to the bath and thermalization effects that, generally, require simultaneous propagation of the bath modes and the system Hamiltonian resulting in the need for extensive computational time.
Another limitation of the method is that the modes treated quantum mechanically should be well defined. OH-stretch vibrations are an example of such well-defined modes as the OH-bond under normal circumstances does not break at ambient temperatures. On the contrary, hydrogen bonds form and break all the time and such low-frequency volatile bond vibrations are not suited for this quantum-classical treatment as their nature is dynamically changed on the short time-scale of the spectroscopy. If a vibration of interest is strongly coupled with many other modes, this also poses a problem as all the strongly coupled modes would need to be included in the system Hamiltonian for an accurate description requiring extensive parameterization and long simulation times. In such case, the methods discussed in Sec. III of this paper may be worth considering.
V. PERSPECTIVES AND A FEW CONCLUDING REMARKS
While 2D IR spectra contain valuable information on couplings, anharmonicity, and frequency correlation, more information can be obtained through higher-order methods generating multidimensional spectra. Three-dimensional infrared spectroscopy has been explored experimentally202 and simulated203,204 with a method based on that described in Sec. IV. Still the interpretation of such spectra remains very challenging. Infrared methods can also be combined with electronic spectroscopy205,206 to reveal the vibrations on the electronic state or to utilize the longer lifetime of the electronic state to probe reaction dynamics. With fluorescence detection,207 2D IR spectroscopy and microscopy would also be possible.208
In the NISE method described in Sec. IV, an essential step is the modeling and parameterization of the multi-chromophore Hamiltonian from electronic structure calculations on their static structures. By contrast, the classical mechanical approach in Sec. III relies on the accurate description of the molecular vibrational property in the context of classical dynamics. This often requires the use of sophisticated molecular mechanics potential models21,209 or quantum mechanical potentials. In this regard, recently proposed ab initio theories of vibrational solvatochromism135–137 and direct QM/MM or full ab initio MD simulations of vibrational spectra210–213 can be fruitfully adapted for either approach to improve their efficiency or accuracy.
2D vibrational spectra can be calculated using the exact hierarchical equations of motion (HEOM) approach.80,214–215 It should be, however, noted that this HEOM approach is based on the assumption of a spectral density description of the bath, which means that the bath is modelled as a collection of independent harmonic oscillators, i.e., the Caldeira-Leggett quantum dissipative bath model.216 Under this assumption, the quantum correlation between the system and the bath can be fully accounted for. However, this method scales rather unfavorably, limiting its application to relatively small systems even though a number of approximations have been developed to improve efficiency.217 Another limitation is that the method requires the parameterization of the spectral density.
In summary, we have presented a brief introduction to the fundamental theory and computational methods used to numerically simulate linear and nonlinear vibrational response functions and the corresponding spectra. One of the widely used methods is to approximate the associated quantum mechanical nonlinear response function with invoking classical approximations so that the linear and nonlinear vibrational response functions can be evaluated by using classical MD simulation trajectories of either equilibrium or non-equilibrium systems. The second method of choice is to solve the time-dependent Schrödinger equation of coupled oscillators where their frequencies, coupling constants, and anharmonicities that are fluctuating in time due to the system-bath interactions can be obtained from independent computational methods. Here, a few representative cases were presented and discussed. Over the past decade, a variety of 2D vibrational spectroscopic techniques that utilize femtosecond IR, THz, and visible pulses have been developed and used to study large-scale delocalized modes in condensed phases and ultrafast reaction dynamics during chemical and biological reactions. Thus, we anticipate that the computational and theoretical methods for accurate calculations of coherent multidimensional vibrational spectroscopic signals of coupled multi-chromophore systems could be incorporated into various MD simulation and quantum chemistry calculation program packages to help the experimentalists in this and related areas to interpret their experimental results and to further understand underlying principles, mechanisms, and functions of materials and molecules in condensed phases.
ACKNOWLEDGMENTS
This work was supported by IBS-R023-D1 (M.C.) and JSPS KAKENHI under Grant No. JP16H02254 (S.S.).
The authors declare no competing financial interests.