HEOM-Based Numerical Framework for Quantum Simulation of Two-Dimensional Vibrational Spectra in Molecular Liquids (HEOM-2DVS)

Abstract

The multi-mode anharmonic Brownian motion model offers a universal framework for simulating molecular vibrations in condensed phases. When vibrational energy surpasses thermal excitation, quantum effects become significant, necessitating a rigorous treatment of system-bath entanglement. The hierarchical equations of motion (HEOM) provide a powerful methodology for simulating such open quantum systems. In this context, two-dimensional vibrational spectroscopy (2DVS) constitutes a powerful probe for elucidating the complex dynamics of molecular processes, both experimentally and theoretically. This work introduces a computational implementation, HEOM-2DVS, for treating non-Markovian open quantum dynamics that encompass energy relaxation, dephasing, thermal excitation, and related processes arising from non-perturbative and nonlinear interactions between selected vibrational modes and their thermal environments. To validate the theoretical framework, we computed

Figures (6)

• Figure 1: Optical Liouville pathways in 2D vibrational spectroscopy for rephasing contribution. In each diagram, the left-hand line depicts the time evolution of the ket state $\lvert {\bf n} \rangle$, while the right-hand line depicts that of the bra state $\langle {\bf n}' \rvert$. The complex-conjugate pathways, obtained by interchanging the left and right states, are not shown.
• Figure 2: Linear absorption (1DIR) spectrum of water calculated for the two‑mode and three‑mode MAB models using HEOM‑2DVS (quantum) and CHFPE-2DVS (classical). The experimental IR dataIRexp2011 are shown as dashed black curves for comparison. Each spectrum is normalized to its maximum peak intensity. The blue solid curves represent the three‑mode classical result, while the green and red solid curves represent the two‑mode and three‑mode quantum results, respectively.
• Figure 3: 2D correlation IR spectra for the stretching and stretching$\rightarrow$bending motions calculated using the two-mode model, which includes (1) the OH stretching mode ($\omega_1 = 3520$ cm$^{-1}$), and (2) the HOH bending mode ($\omega_2 = 1710$ cm$^{-1}$). Spectral intensities were normalized to the maximum amplitude of streching mode. Because the peak intensity of the lower panels is weaker than in the upper panel, the contour interval was tripled for clarity.
• Figure 4: 2D correlation IR spectra for the bending motion for the two-mode case. As the peak intensity was weaker than that in the upper panel of Fig. \ref{['fgr:2DIR2mode_s']}, the contour interval was tripled for emphasis.
• Figure 5: 2D correlation IR spectra for the stretching and stretching$\rightarrow$bending motions calculated using the three-mode model, which includes (1) the OH stretching mode ($\omega_1 = 3570$ cm$^{-1}$), ($1'$) the OH anti-stretching mode ($\omega_{1'} = 3470$ cm$^{-1}$), and (2) the HOH bending mode ($\omega_2 = 1710$ cm$^{-1}$). The mode-mode coupling among the three modes was set to the strong-coupling values listed in Table \ref{['tab:FitAll2']}. Spectral intensities were normalized to the maximum amplitude of strech peak. Because the peak intensity of the lower panels is weaker than in the upper panel, the contour interval was increased by a factor of ten for emphasis.