Feynman Diagram Description of 2D-Raman-THz Spectroscopy Applied to Water

TL;DR

This paper introduces a simplified Feynman diagram approach to interpret 2D-Raman-THz spectra of water, linking experimental signals to hydrogen-bond vibrations with minimal parameters, enhancing understanding of water's microscopic dynamics.

Contribution

It develops a Feynman diagram-based model for 2D-Raman-THz spectroscopy that explains experimental data using a single vibrational mode with anharmonicity and broadening effects.

Findings

Model explains experimental spectra quantitatively

Hydrogen-bond stretching vibration identified at ~170 cm$^{-1}$

Linewidth dominated by quasi-inhomogeneous broadening with 370 fs correlation time

Abstract

2D-Raman-THz spectroscopy of liquid water, which has been presented recently (Proc. Natl. Acad. Sci. USA 110, 20402 (2013)), directly probes the intermolecular degrees of freedom of the hydrogen-bond network. However, being a relatively new technique, its information content is not fully explored as to date. While the spectroscopic signal can be simulated based on molecular dynamics simulation in connection with a water force field, it is difficult to relate spectroscopic signatures to the underlying microscopic features of the force field. Here, a completely different approach is taken that starts from an as simple as possible model, i.e., a single vibrational mode with electrical and mechanical anharmonicity augmented with homogeneous and inhomogeneous broadening. An intuitive Feynman diagram picture is developed for all possible pulse sequences of hybrid 2D-Raman-THz spectroscopy. It is shown that the model can explain the experimental data essentially quantitatively with a very small set of parameters, and it is tentatively concluded that the experimental signal originates from the hydrogen-bond stretching vibration around 170 cm$^{-1}$. Furthermore, the echo observed in the experimental data can be quantified by fitting the model. A dominant fraction of its linewidth is attributed to quasi-inhomogeneous broadening in the slow-modulation limit with a correlation time of 370 fs, reflecting the lifetime of the hydrogen-bond networks giving rise the absorption band.

Figures

Tables

Table 1. Table 1: Intensities of all non-zero peaks appearing in Eq. ( 11 ).

	Peak
	(a)	(b)	(c)	(d)
Response	($-\omega ,\omega$)	($\omega ,\omega$)	($\omega ,0$)	($\omega ,2\omega$)
$R_{\mathrm{RTT}}^{\alpha }$	0	0	0	0
$R_{\mathrm{TRT}}^{\alpha }$	1	-1	0	0
$R_{\mathrm{TTR}}^{\alpha }$	0	0	1	-1
$R_{\mathrm{RTT}}^{\mu }$	1	-1	1	-1
$R_{\mathrm{TRT}}^{\mu }$	0	0	1	-1
$R_{\mathrm{TTR}}^{\mu }$	1	-1	0	0
$R_{\mathrm{RTT}}^M$	1	3	-3	-1
$R_{\mathrm{TRT}}^M$	1	3	-3	-1
$R_{\mathrm{TTR}}^M$	1	3	-3	-1

Table 2. Table 2: Parameters obtained from the fit of the hydrogen-bond bending vibration and hydrogen-bond stretching vibration.

	stretching	bending
$\omega$ [cm-1]	16511footnotemark: 1	4511footnotemark: 1
$\Delta \omega$ [cm-1]	75	40
${\tau }_c$ [fs]	370	300222As a result of Eq. 21, the peak position of the band does not coincide exactly with $\omega$
$T_1$ [fs]	250	300222As a result of Eq. 21, the peak position of the band does not coincide exactly with $\omega$
${\sigma }_{\mu }/{\sigma }_M$	1.4	0.96
${\sigma }_{\alpha }/{\sigma }_M$	-0.1	0.01

Equations

$$R_{\mathrm{RTT}}$$

$$R_{\mathrm{TRT}}$$

$$R_{\mathrm{TTR}}$$

$${\hat{H}}^{(0)}=\frac{\hslash \omega }{2}\left({\hat{p}}^2+{\hat{q}}^2\right).$$

$${\left({\mathbf{q}}_{\mathbf{H}}\right)}_{ij}\equiv ⟨{\phi }_i|\hat{q}|{\phi }_j⟩=\frac{1}{\sqrt{2}}\left[\sqrt{j}{\delta }_{i+1,j}+\sqrt{i}{\delta }_{i-1,j}\right]$$

$$\mathbf{\alpha }\propto {\mathbf{q}}_{\mathbf{H}}+...$$

$$\mathbf{\mu }\propto {\mathbf{q}}_{\mathbf{H}}+...$$

$$\mathbf{\alpha }\propto {\mathbf{q}}_{\mathbf{H}}+{\sigma }_{\alpha }{{\mathbf{q}}_{\mathbf{H}}}^2+...$$

$$\mathbf{\mu }\propto {\mathbf{q}}_{\mathbf{H}}+{\sigma }_{\mu }{{\mathbf{q}}_{\mathbf{H}}}^2+...$$

$$\displaystyle(\boldsymbol{q_{H}}^{2})_{ij}=\tfrac{1}{2}\Big{[}(2i+1)\delta_{i,j}$$

$$\displaystyle\sqrt{(j+1)(j+2)}\delta_{i,j+2}\Big{]},$$

$$\hat{H}={\hat{H}}^{(0)}+{\sigma }_{\mathrm{M}}\hslash \omega q^3$$

$$\displaystyle\left(\boldsymbol{q_{M}}\right)_{ij}=\tfrac{1}{2}\Big{[}-3(2i+1)\delta_{i,j}$$

$$\displaystyle\sqrt{(j+1)(j+2)}\delta_{i,j+2}\Big{]},$$

$$\mathbf{\alpha }$$

$$\approx {\mathbf{q}}_{\mathbf{H}}+{\sigma }_{\mathrm{M}}{\mathbf{q}}_{\mathbf{M}}+{\sigma }_{\alpha }{{\mathbf{q}}_{\mathbf{H}}}^2$$

$$\mathbf{\mu }$$

$$\approx {\mathbf{q}}_{\mathbf{H}}+{\sigma }_{\mathrm{M}}{\mathbf{q}}_{\mathbf{M}}+{\sigma }_{\mu }{{\mathbf{q}}_{\mathbf{H}}}^2$$

$$R_{\mathrm{RTT}}\propto {\sigma }_{\alpha }{R}_{\mathrm{RTT}}^{\alpha }+{\sigma }_{\mu }{R}_{\mathrm{RTT}}^{\mu }+{\sigma }_{\mathrm{M}}{R}_{\mathrm{RTT}}^{\mathrm{M}}$$

$$R_{\mathrm{TRT}}\propto {\sigma }_{\alpha }{R}_{\mathrm{TRT}}^{\alpha }+{\sigma }_{\mu }{R}_{\mathrm{TRT}}^{\mu }+{\sigma }_{\mathrm{M}}{R}_{\mathrm{TRT}}^{\mathrm{M}}$$

$$R_{\mathrm{TTR}}\propto {\sigma }_{\alpha }{R}_{\mathrm{TTR}}^{\alpha }+{\sigma }_{\mu }{R}_{\mathrm{TTR}}^{\mu }+{\sigma }_{\mathrm{M}}{R}_{\mathrm{TTR}}^{\mathrm{M}}$$

$$R_{\mathrm{RTT}}^{\alpha }=$$

$$R_{\mathrm{TRT}}^{\alpha }=$$

$$R_{\mathrm{TTR}}^{\alpha }=$$

$$R_{\mathrm{RTT}}^{\mu }=$$

$$-\mathrm{Tr}\left\{{{\mathbf{q}}_{\mathbf{H}}}^2(t_1+t_2)\left[{\mathbf{q}}_{\mathbf{H}}(t_1),\left[{\mathbf{q}}_{\mathbf{H}}(0),{\rho }_{\mathrm{eq}}\right]\right]\right\}$$

$$R_{\mathrm{TRT}}^{\mu }=$$

$$-\mathrm{Tr}\left\{{{\mathbf{q}}_{\mathbf{H}}}^2(t_1+t_2)\left[{\mathbf{q}}_{\mathbf{H}}(t_1),\left[{\mathbf{q}}_{\mathbf{H}}(0),{\rho }_{\mathrm{eq}}\right]\right]\right\}$$

$$R_{\mathrm{TTR}}^{\mu }=$$

$$-\mathrm{Tr}\left\{{\mathbf{q}}_{\mathbf{H}}(t_1+t_2)\left[{{\mathbf{q}}_{\mathbf{H}}}^2(t_1),\left[{\mathbf{q}}_{\mathbf{H}}(0),{\rho }_{\mathrm{eq}}\right]\right]\right\}$$

$$R_{\mathrm{RTT}}^{\mathrm{M}}=$$

$$-\mathrm{Tr}\left\{{\mathbf{q}}_{\mathbf{H}}(t_1+t_2)\left[{\mathbf{q}}_{\mathbf{M}}(t_1),\left[{\mathbf{q}}_{\mathbf{H}}(0),{\rho }_{\mathrm{eq}}\right]\right]\right\}$$

$$-\mathrm{Tr}\left\{{\mathbf{q}}_{\mathbf{H}}(t_1+t_2)\left[{\mathbf{q}}_{\mathbf{H}}(t_1),\left[{\mathbf{q}}_{\mathbf{M}}(0),{\rho }_{\mathrm{eq}}\right]\right]\right\}$$

$$=R_{\mathrm{TRT}}^{\mathrm{M}}=R_{\mathrm{TTR}}^{\mathrm{M}},$$

$$\displaystyle R_{\mathrm{RTT}}(t_{1},t_{2})\propto\Theta(t_{1},t_{2})\big{(}$$

$$+$$

$$\displaystyle R_{\mathrm{TRT}}(t_{1},t_{2})\propto\Theta(t_{1},t_{2})\big{(}$$

$$+$$

$$\displaystyle R_{\mathrm{TTR}}(t_{1},t_{2})\propto\Theta(t_{1},t_{2})\big{(}$$

$$+$$

$${\Gamma }_{-\omega ,\omega }(t_1,t_2)=$$

$${\Gamma }_{\omega ,\omega }(t_1,t_2)=$$

$${\Gamma }_{\omega ,0}(t_1,t_2)=$$

$${\Gamma }_{\omega ,2\omega }(t_1,t_2)=$$

$$T_1^{(n,n)}$$

$$T_2^{(n,m)}$$

Full text

I Introduction

THz spectroscopy of liquids probes their intermolecular modes directly. Applied to liquid water, it reveals three distinct bands, which are assigned to hydrogen bond bending modes (50 cm*-1*), hydrogen bond vibrations (170 cm*-1*), and librations, which are hindered rotations of water molecules (600 cm*-1*).Bertie1996 These spectra can give us a picture of the motions in water, however, vibrational coherences are very short lived in water due to the fast, chaotic dynamics of the hydrogen bonding network, causing very broad and blurred bands. This limits the amount of information that can effectively be extracted from these spectra.

Extending the spectroscopy into two dimensions can thin out the information, and thus increase the amount of accessible information, and can furthermore disentangle homogeneous from inhomogeneous broadening. tan93 ; okumura97a There has been significant effort to extend the very concept of 2D spectroscopy into the THz regime, either in the form of 2D Raman spectroscopy VBout1991 ; Inaba1993 ; Muller1993 ; Tokmakoff1997 ; Blank1999 ; Blank2000 ; Kaufman2002 ; Golonzka2000 ; Kubarych2003 ; Li2008 , 2D-THz spectroscopy,Kuehn2009 ; Kuehn2011 ; Kuehn2011a ; Elsaesser2015 ; Lu2016 ; Somma2016 or hybrid 2D-Raman-THz methods.Finneran2016 ; Finneran2017 ; Savolainen2013 ; Shalit2017 Due to technical limitations, the latter is as of now the only 2D spectroscopy in the THz range that has been successfully applied to water and aqueous salt solutions.Savolainen2013 ; Shalit2017 ; Berger2018 Of particular interest in these experiments is the observation of a very short-lived echo, whose decay time has been related to the relevant time-scales in water. However, a detailed understanding of the 2D-Raman-THz response is still lacking.

Unlike conventional 2D IR spectroscopy,hamm11 2D-Raman spectroscopy as well as 2D-Raman-THz spectroscopy are described by second-order perturbation theory, despite the fact that it is a 3rd-order response with regard to the electrical fields. This is since the Raman process is electronically non-resonant and thus instantaneous, and it is assumed that the two field interactions giving rise to the Raman process occur at the same time. The three instead of four interactions of the system (including the emission process) bring about that one has to induce a two-quantum transition, which would be forbidden in the harmonic approximation. This forbidden transition is a bottleneck of the signal, which is why its very cause (i.e., electrical or mechanical anharmonicity) determines the shape of the signal. Thus, the two-quantum transition must be an integral part of any model used to describe 2D-Raman-THz spectroscopy.

Time dependent perturbation theory is the starting point to calculate the 2D-Raman-THz response. The response functions for the three different time-orderings are:Cho1999

[TABLE]

where $\hat{\mu}$ is the dipole operator, $\hat{\alpha}$ the polarizability operator, $t_{1}$ the time between the first and the second laser pulse interacting with the sample, and $t_{2}$ the time from the second laser pulse to the detection process. We have concentrated on the Raman-THz-THz (RTT) and the THz-Raman-THz (TRT) pulse sequences,Hamm2012 ; Hamm2012a ; Savolainen2013 ; Hamm2014 ; Hamm2017 ; Shalit2017 ; Berger2018 where the detection step measures an emitted THz field, while Blake and coworkers looked at the THz-THz-Raman (TTR) pulse sequence with a Raman process for detection.Finneran2016 ; Finneran2017

The theory of 2D-Raman-THz spectroscopy is similar to 2D-Raman spectroscopy,tan93 ; Steffen1996 ; Steffen1998 ; saito98 ; okumura97a ; ma00 ; Jansen2000 ; Okumura2003 ; saito06 ; Hasegawa2006 as well as hybrid IR-Raman techniques.Cho1999 ; Cho2000 ; Zhao2000 ; Guo2009 ; Grechko2018 In addition, a fair share of theory has been published tailored specifically for 2D-Raman-THz spectroscopy.Hamm2012 ; Hamm2012a ; Hamm2014 ; Ito2014 ; Ito2015 ; Ikeda2015 ; Pan2015 ; Ito2016 ; Ito2016b Since typical THz experiments work in a frequency range equivalent to $k_{B}T$, the response can be derived in the classical limit from molecular dynamics (MD) simulation.Hasegawa2006 This approach appears to be the method of choice for complicated systems like water, since basically all effects, apart from possible quantum effects,Berger2018 are captured implicitly by a MD force field, including anharmonicities, mode coupling, chemical exchange, and orientational averaging. It has been shown that the MD approach reveals responses, which strongly depend on the force field used, especially on the description of polarizability, albeit in a rather nonintuitive and indirect way.Hamm2014 ; Ito2015 These MD simulations are largely “black-box” and it is difficult to disentangle the contributions to the 2D-Raman-THz response and to relate spectroscopic signatures to the underlying microscopic features of a water force field.

In order to learn more about the microscopic mechanism giving rise to the 2D-Raman-THz response, we take a completely different approach here and start from an as simple as possible model, i.e., a single vibrational mode augmented with homogeneous and inhomogeneous broadening. To that end, we follow the conceptual framework introduced in Ref. Okumura2003, for the description of 2D Raman spectroscopy, starting from a quantum-mechanical harmonic oscillator in an eigenstate representation and adding electrical and/or mechanical anharmonicity. The response functions Eq. 1 can then be depicted in a very intuitive way in terms of Feynman diagrams, from which one can directly read off the peak position in a 2D spectrum. We will see that the model can explain the experimental data essentially quantitatively with a small set of parameters.

II Model System

II.1 Zero-Order Description: Harmonic oscillator

Due to its simplicity and the fact that it is often a very good approximation for molecular vibrations, the harmonic oscillator is an obvious starting point to set the stage:

[TABLE]

Analytical solutions of energy levels and eigenfunctions exist. From the thermal population of eigenstates, the equilibrium density matrix $\rho_{\mathrm{eq}}$ in Eq. 1 can readily be constructed. Furthermore, the dimensionless position operator $\hat{q}=\sqrt{\frac{m\omega}{\hbar}}\hat{x}$ can be expressed in an harmonic oscillator eigenstate basis $|\varphi_{i}\rangle$

[TABLE]

(where the subscript in $\boldsymbol{q_{H}}$ stands for “harmonic”), which can be seen from the common creation ($\hat{b}$) and annihilation ($\hat{b}^{\dagger}$) operator formalism with $\hat{q}=1/\sqrt{2}(\hat{b}^{\dagger}+\hat{b})$.Atkins2010 With that, the dipole and polarizability operators can be expanded in $\boldsymbol{q_{H}}$:

[TABLE]

The proportionality factors are irrelevant for the purpose of this paper, as they give rise to an overall intensity of the 2D-Raman-THz signal, which however is typically not determined experimentally in an absolute sense. The $\delta_{i+1,j}$ and $\delta_{i-1,j}$ terms in $\boldsymbol{q_{H}}$ couple states $i$ and $i\pm 1$, which leads to the well-known $i\rightarrow i\pm 1$ selection rules of the harmonic oscillator. This level of description is often sufficient to describe linear (1D) THz or Raman spectroscopy.

The $i\rightarrow i\pm 1$ selection rules cause coherence pathways, in which the density matrix is alternating between population and coherence states. As one starts from a thermal population state $\rho_{\mathrm{eq}}$, an even number of interactions would be needed to return to a population state after the emission of a signal. Therefore, due to the odd number of interactions in 2D-Raman-THz spectroscopy, the harmonic oscillator together with Eq. 4 would predict a vanishing signal. The appearance of a signal requires zero- or two-quantum transitions, which can be accomplished by breaking the symmetry of the system; either by considering nonlinearity of $\boldsymbol{\mu}$ and $\boldsymbol{\alpha}$ (electrical anharmonicity), or by perturbing the potential of the oscillator (mechanical anharmonicity). Softening the selection rules results in a bottleneck for the response, hence the signal shape is very sensitive to the specific nature of the anharmonicity, which thus is a crucial aspect of the modelling of 2D-Raman-THz spectra.

II.2 Electrical Anharmonicity

Electrical anharmonity is introduced by considering higher order terms in equation Eq. (4): Okumura2003

[TABLE]

with dimensionless smallness parameters $|\sigma_{\mu}|\ll 1$ and/or $|\sigma_{\alpha}|\ll 1$. In an harmonic eigenstate basis, the quadratic term is:

[TABLE]

which again can be seen from the creation/annihilation operator formalism.Atkins2010 In analogy to an electrical quadrupole transition, the $\delta_{i,j}$-term allows for zero-quantum transitions $i\rightarrow i$, and the $\delta_{i,j\pm 2}$-terms for two-quantum transitions $i\rightarrow i\pm 2$.

II.3 Mechanical Anharmonicity

Mechanical anharmonicity (which has not been considered in Ref. Okumura2003, ) breaks the symmetry by modifying the potential energy function. We consider a cubic anharmonic oscillator by adding $\sigma_{\mathrm{M}}\hbar\omega q^{3}$ to the harmonic Hamiltonian:

[TABLE]

with dimensionless smallness parameter $|\sigma_{\mathrm{M}}|\ll 1$. The eigenstates of the anharmonic oscillator $|\varphi_{i}^{(anh)}\rangle$ can be calculated perturbatively.Claude1991 Although the position operator $\hat{q}$ per se is not affected by this addition, its matrix representation in an anharmonic eigenstate basis is. It can be shown (see Appendix A) that this matrix representation takes the form $\langle\varphi_{i}^{(anh)}|\hat{q}|\varphi_{j}^{(anh)}\rangle=(\boldsymbol{q_{H}})_{ij}+\sigma_{\mathrm{M}}(\boldsymbol{q_{M}})_{ij}$ with:

[TABLE]

where the subscript “M” stands for “mechanic anharmonicity”. $\boldsymbol{q_{M}}$ has the same structure as $\boldsymbol{q_{H}}^{2}$ (Eq. 6) with the $\delta_{i,j}$-term causing zero-quantum transitions and the $\delta_{i,j\pm 2}$-terms causing two-quantum transitions, however the prefactors determining the transition probabilities and the signs of the peaks are different.

For a system with mechanical and electrical anharmonicity at the same time, the operators are:

[TABLE]

where we neglected in the second step terms that are higher than first order in any of the smallness parameters $\sigma_{\mu}$, $\sigma_{\alpha}$ and $\sigma_{\mathrm{M}}$. For the same reason, we also neglect coherence pathways that contain more than one forbidden transition.

II.4 Feynman Diagrams

Based on Eq. 9, the response functions of Eq. 1 can be separated into:

[TABLE]

with

[TABLE]

which allows us to discuss the contributions from electrical ($\sigma_{\alpha}$,$\sigma_{\mu}$) and mechanical anharmonicity ($\sigma_{\mathrm{M}}$) separately. To this end, we expand the nested commutators in Eqs. 11 and represent all terms in the form of Feynman diagrams (Figs. 1, 2 and 3). The corresponding 2D peaks are also summarized in Tab. 1 together with their intensities, which originate from products of the prefactors of the corresponding $\delta$-terms in Eqs. 3, 6 and 8, and their sums when contributions overlap in frequency space. The frequency positions of the in total four peaks, ($-\omega,\omega$), ($\omega,\omega$), ($\omega,0$), and ($\omega,2\omega$), are the same as in Ref. Okumura2003, , but their intensities differ since the selection rules of the three pulse sequences of 2D-Raman-THz spectroscopy differ from those in 2D-Raman spectroscopy. The various spectra are discussed in the following.

In $R_{\mathrm{RTT}}^{\alpha}$, $R_{\mathrm{TRT}}^{\alpha}$, and $R_{\mathrm{TTR}}^{\alpha}$ (Fig. 1), the nonlinearity of $\boldsymbol{\alpha}$ allows the Raman interaction to induce a zero- or a two-quantum transition via the contribution of $\boldsymbol{q_{H}}^{2}$. Starting with $R_{\mathrm{RTT}}^{\alpha}$, four Feynman diagrams can be constructed (Fig. 1, top left). In this case the first Raman interaction induces a forbidden transition, i.e. a zero-quantum or a two-quantum transition. The pathways contributing to peak (e) start with a zero-quantum transition from either the left and the right, which however lead to two equivalent contributions of opposite sign that perfectly cancel out. Peak (f), on the other hand, starts with a two-quantum transition from the left. From there, the system goes into a $|1\rangle\langle 0|$ coherence by the second interaction from the left, or into $|2\rangle\langle 1|$ coherence via an interaction from the right. For a harmonic potential, the energy levels are equally spaced and both coherences oscillate with the same frequency. If we furthermore assume that the 0-1 and the 1-2 dephasing times are the same (see Eq. 14 below),Okumura2003 both contributions again cancel, and we have $R_{\mathrm{RTT}}^{\alpha}=0$ altogether.

In contrast, two peaks appear in the TRT pulse sequence. For peak (a), the first interaction brings the system into a $|0\rangle\langle 1|$ coherence, and the second Raman interaction does a two-quantum transition to bring the system into the $|2\rangle\langle 1|$ coherence, hence the peak appears at $(-\omega,\omega)$. Peak (b), in contrast, is diagonal at $(\omega,\omega)$, since the second interaction does a zero-quantum transition, which does not affect the state. The zero-quantum transition can either act from the left or the right side, giving rise to two contributions with opposite sign. However, these contributions do not fully cancel, since the interaction from left is weighted by the $(\boldsymbol{q_{H}}^{2})_{1,1}=-3/2$ element, and the interaction from right by $(\boldsymbol{q_{H}}^{2})_{0,0}=+1/2$. Finally, the TTR sequence has its forbidden transition in the emission step with zero-quantum coherences for peak (c), and two-quantum coherences for peak (d). Peak (c), with $\omega_{2}=0$, appears together with its conjugate complex (c∗). The conjugate complex of all other peaks show up at negative frequencies $\omega_{2}$, but the negative $\omega_{2}$-half-space is not shown in Fig. 1 for clarity.

There are more peaks in $R_{\mathrm{RTT}}^{\mu}$, $R_{\mathrm{TRT}}^{\mu}$, and $R_{\mathrm{TTR}}^{\mu}$ (Fig. 2), since there are two interactions with $\boldsymbol{\mu}$ and the forbidden transition can occur with either one of it. However, Eq. 11 shows that each term of $R^{\mu}$ can be written as a sum of two terms that we already discussed in the context of $R^{\alpha}$, e.g. $R_{\mathrm{RTT}}^{\mu}=R_{\mathrm{TRT}}^{\alpha}+R_{\mathrm{TTR}}^{\alpha}$. The corresponding 2D spectra in Fig. 2 therefore are overlays of two spectra each from Fig. 1.

Finally, the responses $R_{\mathrm{RTT}}^{\mathrm{M}}$, $R_{\mathrm{TRT}}^{\mathrm{M}}$ and $R_{\mathrm{TTR}}^{\mathrm{M}}$ originating from mechanical anharmonicity are shown in Fig. 3. The term $\sigma_{\mathrm{M}}\boldsymbol{q_{M}}$ contributes to $\boldsymbol{\mu}$ and $\boldsymbol{\alpha}$ equally, and the forbidden transition can occur at any of the tree interactions. Since $\boldsymbol{\mu}=\boldsymbol{\alpha}$ for this case, the order of interactions does not matter, and $R_{\mathrm{RTT}}^{\mathrm{M}}=R_{\mathrm{TRT}}^{\mathrm{M}}=R_{\mathrm{TTR}}^{\mathrm{M}}$. The corresponding 2D spectrum is in essence an overlay of all three spectra shown in Fig. 1, except for the fact that the features that include zero-quantum transitions (b,c) change signs and amplitudes due to the different diagonal element of $\boldsymbol{q_{M}}$. Furthermore, peak (f) is now nonzero due to the anharmonic shift of the energy levels.

As in Ref. Okumura2003, , we considered so far only Feynman diagrams starting from the $|0\rangle\langle 0|$ element of the density matrix in Figs. 1, 2 and 3, implicitly assuming a temperature $T$=0 K. It can however be shown that the overall response function does not depend on the starting level, despite the fact that more Feynman diagrams come into play, and hence is in fact temperature independent (see Appendix B for details). For that we however need to assume that the lineshape functions depend on level differences only (see Eq. 14 below), and that the effect of mechanical anharmonicity is on the transitions probabilities only (via $\boldsymbol{q_{M}}$), while we keep the energy spectrum equidistant. Fig. 3 shows that this is a good approximation. That is, the anharmonic shift $\omega_{21}-\omega_{10}$ is small, and for any reasonable, not too narrow spectral linewidth both contributions giving rise to peak (f) largely cancel each other. In comparison, the effect of mechanical anharmonicity on the transition probabilities, giving rise to peaks (a)-(e), is much larger. Within that approximation, quantum and classical response functions are in fact identical, and simple analytic expressions can be derived for the three response functions by collecting all terms (without that approximation, an explicit sum over Boltzmann-weighted initial states would be needed that converges extremely slowly for $\hbar\omega\ll k_{B}T$):

[TABLE]

where $\Theta(t_{1},t_{2})$ is the Heaviside function ensuring that the response vanishes for $t_{1}<0$ or $t_{2}<0$. In second order perturbation theory, the complex conjugate of a term can be generated by two permutations in the commutator. Therefore, the two complex conjugate components have the same sign and reveal $\cos$-functions upon addition.

We augment these response functions with dephasing terms $\Gamma(t_{1},t_{2})$ in essence along the lines of Ref. Steffen1998, , but replacing the part originating from vibrational $T_{1}$ relaxation with Eq. 36, rather than Eq. 37, from Ref. Okumura2003, . That is, each response function in Eq. 12 contains one term each for peaks at $(-\omega,\omega)$, $(\omega,\omega)$, $(\omega,0)$ and $(\omega,2\omega)$, respectively, which are multiplied with the following damping functions:

[TABLE]

As discussed in Ref. Okumura2003, , the damping terms caused by vibrational relaxation depend on quantum numbers according to:

[TABLE]

The lineshape function $g(t)$ is:

[TABLE]

whose origin are Gaussian fluctuations of the transition frequency $\omega(t)$ with standard deviation $\Delta\omega$ and correlation time $\tau_{c}$:

[TABLE]

where $\delta\omega(t)=\omega(t)-\langle\omega(t)\rangle$. In the limes $\Delta\omega\tau_{c}\ll 1$, Eq. 15 causes pure dephasing with $T^{*}_{2}=(\Delta\omega^{2}\tau_{c})^{-1}$, in the limes $\Delta\omega\tau_{c}\gg 1$ inhomogeneous dephasing with a Gaussian lineshape with width $\Delta\omega$.

II.5 Inhomogeneous Broadening and Echoes

In any of the responses, peak (a) with frequencies of opposite signs $(-\omega,\omega)$ is a rephasing pathway,Hamm2017 as a coherence that dephased during the first time period might rephase during the second, provided the band is inhomogeneously broadened and maintains some amount of memory on the oscillation frequency. Rephasing requires an “inversion of coherence”, and the only possibility to achieve that are coherence pathway that start with a $|0\rangle\langle 1|$ coherence and are then brought into a $|2\rangle\langle 1|$ coherence by the second pulse (see Figs. 1, 2 and 3). In a time-domain representation, rephasing will generate an echo along $t_{1}=t_{2}$; in a frequency-domain representation tilted 2D lineshapes. Peak (b), in contrast, is a non-rephasing pathway.

To explore the appearance of echoes (Fig. 4), we choose lineshape parameters for the hydrogen-bond stretching vibration as listed in Tab. 2 (see below), except of unrealistically long values for $\tau_{c}=$100 ps and $T_{1}=1~{}$ps. Echoes along $t_{1}=t_{2}$ are seen for $R_{TRT}^{\alpha}$, $R_{RTT}^{\mu}$, and $R_{TTR}^{\mu}$, i.e., coherence pathways with their two-quantum transition for the second pulse to induce the required two-quantum transition. Experimentally, one can distinguish between $R_{RTT}$, $R_{TRT}$, and $R_{TTR}$ by the choice of a pulse sequence, hence when observing an echo, one gets an handle on determining the major source of anharmonicity. The vertical features in $R_{TTR}^{\alpha}$, $R_{RTT}^{\mu}$ and $R_{TRT}^{\mu}$ originate from peak (c), which does not experience any inhomogeneous dephasing in the $t_{2}$ direction, since its frequency in $\omega_{2}$ is always zero.

III Comparison to Water Experiments

III.1 Instrument Response Function

We now explore to what extent the model derived here can explain the experimental data of liquid water from Ref. Savolainen2013, . We start with noting that this experiment measured the $R_{\mathrm{RTT}}$ and $R_{\mathrm{TRT}}$ responses by interchanging the timing between the Raman and THz-pump pulses, but not the $R_{\mathrm{TTR}}$ response, which would require a different detection scheme.Finneran2016 ; Finneran2017 With the particular arrangement of the delay lines in the experiment of Ref. Savolainen2013, , we measured:

[TABLE]

where time $t_{1}$ is from the Raman pump pulse to the THz pump pulse and time $t_{2}$ from the THz pump pulse to detection (which is why $t_{2}-t_{1}$ appears as argument of $R_{\mathrm{TRT}}$ where the THz pump comes before the Raman pump, while $t_{2}$ is the time from the second Raman interaction to detection in the definition of $R_{\mathrm{TRT}}$ in Eq. 12). Furthermore, the experimental signal is obtained from the response functions by a convolution with the instrument response function (IRF):

[TABLE]

which is calculated from the THz field $E_{\mathrm{THz}}$ and the envelope of the Raman pulse $I_{\mathrm{Raman}}$:

[TABLE]

Eq. 19 is an idealized expression to illustrate the basic idea; in the real experiment the IRF contains in addition a transfer function describing how the emitted THz field reshapes as it propagates from the sample to the detection crystal, as well as a correction for a non-perfect Gouy phase (see Refs. Savolainen2013, ; Ahmed2014, for details).

It is illustrative to look at the IRF in the frequency domain (Fig. 5), as it shows which one of the peaks in the 2D responses (Figs. 1, 2 and 3) the experiment is sensitive to. The IRF spectrum for the RTT pulse sequence (Fig. 5a) has two nodal lines, one at $\omega_{2}=0$, which is caused by the derivative $d/dt_{2}$ in Eq. 19, and one at $\omega_{1}=\omega_{2}$, which caused by the fact that the THz pump pulse does not have any DC component (i.e., $\int E_{\mathrm{THz}}(t)dt=0$). Peaks (b) and (c) lie on these nodal lines, as they would contain zero-quantum transitions for either the THz pump pulse or the THz emission process. We will nevertheless see these peaks to a certain extent since the damping in water is very fast, hence they have a significant spectral width that extends into regions where the IRF is non-zero. Furthermore, a smaller bandwidth is observed in the rephasing quadrant ($\omega_{1}<0$, $\omega_{2}>0$) than in the non-rephasing quadrant ($\omega_{1},\omega_{2}>0$), reflecting the fact that the peak at $\omega_{2}=-\omega_{1}$ includes a two-quantum transition which requires a higher bandwidth of the THz pulse.

For the $R_{\mathrm{TRT}}$ response, the IRF has to be transformed to account for the time shift applied in Eq. 17 (Fig. 5b). The effect of a time transform $t_{1}=-t_{1}^{\prime}$ and $t_{2}=t_{2}^{\prime}-t_{1}^{\prime}$ on the frequency domain can be derived by inspecting the effect on a single basis function of the Fourier transformation:

[TABLE]

with $\omega_{1}^{\prime}=-\omega_{1}-\omega_{2}$ and $\omega_{2}^{\prime}=\omega_{2}$. The transformed instrument response has nodal lines at $\omega_{1}=0$ and $\omega_{2}=0$. Hence, in this case, the diagonal peak (b) lies in a region where the IRF is large, since the envelope of the Raman pulse has a zero-frequency component.

III.2 1D Spectra and 2D Response

Fig. 6 shows the 1D Raman spectrum of water in the relevant frequency range, which contains two spectroscopic features: the hydrogen-bond bending vibration around 50 cm*-1* and the hydrogen-bond stretch vibration around 170 cm*-1* (the librational mode around 600 cm*-1* is not considered here since it is completely ouside our experimental observation window). Also shown in Fig. 6 are fits to these two bands, assuming a linear response function

[TABLE]

The parameters are listed in Tab. 2 and will be justified later based on the fit of the echo observed in the 2D response.

Starting with the hydrogen-bond stretching vibration, Fig. 7 (left column) shows the corresponding 2D results in the time domain, i.e., the convolution of Eq. 17 with the IRF. The lineshape parameters of Tab. 2 were assumed, except for $\tau_{c}=$100 ps and $T_{1}=1~{}$ps, which in a first step was chosen unrealistically slow in order to more clearly demonstrate the resulting echoes (i.e., the same parameters as in Fig. 4). Each one of the three types of anharmonicities give rise to a different 2D response. The echoes shown in Fig. 4 survive the convolution and are marked by arrows in Fig. 7. In the case of a nonlinear dipole moment, the rephasing peak is in $R_{\mathrm{RTT}}^{\mu}$ and the echo is visible in the (upper right) RTT quadrant along $t_{1}=t_{2}$. On the other hand, for $S_{\alpha}$ with the nonlinearity in the polarizability, the rephasing peak is in the $R_{\mathrm{TRT}}^{\alpha}$ response and the echo appears in the upper left quadrant along the $t_{2}=-2t_{1}$ due to the time transformation in Eq. 17. Finally, $S_{\mathrm{M}}$ is a mix of both (albeit with different signs), since mechanical anharmonicity allows for all coherence pathways at the same time. Fig. 7 (right column) shows the same, but now for $\tau_{c}=370$ fs and $T_{1}$=250 fs (which will be justified later based on a fit to the experimental data).The echoes are now masked by the instrument response function.

Modelling the hydrogen-bond bending vibration around 50 cm*-1* on the same footing reveals similar results, see Fig. 8. The echoes are clearer in this case, even when using realistic parameters for $\tau_{c}$ and $T_{1}$ (Fig. 8, right column). In addition the signal intensity is significantly larger (Fig. 7 has been upscaled by a factor 2), reflecting the the limited bandwidth of the THz pulses in the experiment, that peak at 50 cm*-1* but only partially cover the 170 cm*-1* band.Savolainen2013

III.3 Fit of the Water Response

The free parameters of the model are $\sigma_{\mu}$, $\sigma_{\alpha}$, $\sigma_{\mathrm{M}}$, $T_{1}$, $\tau_{c}$, $\omega$ and $\Delta\omega$. In an iterative process, we varied the dephasing parameter $T_{1}$, $\tau_{c}$ to reproduce the 2D-response, and $\omega$ as well as $\Delta\omega$ to reproduce the position and width of the hydrogen-bond bending vibration in the 1D spectrum (Fig. 6). In addition, small variations of the delay-time zeros and the correction for the Gouy phase were allowed, which are not very accurately defined in the experiment.Savolainen2013 Once these nonlinear parameters are fixed, the 2D responses of Figs. 7 and 8 (right column) can be considered a basis in which the experimental response is expanded in order to minimize the RMSD between experimental and fitted spectrum. Fig. 9 shows that this procedure results in a remarkably good agreement with the experimental data, despite the simplicity of the model and the small number of parameters. Tab. 2 summarizes the resulting parameters for the two modes.

The RMSD of the fit for the hydrogen-bond stretching vibration is smaller by a significant factor 0.7 as compared to that of the hydrogen-bond bending vibration. It has in fact not been possible to find a minimum in the RMSD for the hydrogen-bond bending vibration (see footnote in Tab. 2), and in addition, the shift in $t_{2}$ required to obtain the best fit for hydrogen-bond bend vibration is larger (i.e., $\Delta t_{2}$=170 fs) than what we would think is its experimental uncertainty (the corresponding value for the hydrogen-bond stretching vibration is $\Delta t_{2}=35$ fs). We therefore suggest that the observed experimental signal originates predominantly from the hydrogen-bond stretching vibration, despite the fact its response is reduced to a certain extent as it is at the very edge of the experimentally accessible frequency window. We have not considered scenarios in which both bands contribute in parallel, or even more so, couple with each other, because the number of fitting parameters would be too large and the fitting problem would be under-determined.

Fig. 10 plots the RMSD of the fit for the hydrogen-bond stretching vibration as a function of $\tau_{c}$ and $T_{1}$, revealing that both parameters are not strongly correlated. As these two parameters determine the relative contribution of inhomogeneous vs homogeneous dephasing, we conclude that the signatures of the echo are still present in the 2D response. In that regard it is fortunate that the rephasing peak (a) in $R_{\mathrm{RTT}}^{\mu}$ (Fig. 2, bottom left) and $R_{\mathrm{RTT}}^{\mathrm{M}}$ (Fig. 3, bottom) have the same sign, while all other peaks have opposite sign (and the contribution of $R_{\mathrm{RTT}}^{\mu}$ is negligible, see Tab. 2). Together with the corresponding weights of the two contributions ($\sigma_{\mu}/\sigma_{\mathrm{M}}$=1.4), the rephasing peak will actually dominate in the overall response. That is illustrated in Fig. 11, which shows the response in the the frequency domain and without the convolution with the instrument response function for the same parameters as in in Fig. 9b; $R_{\mathrm{RTT}}$ is dominated by one rephasing peak that is strongly elongated along the diagonal due to the inhomogeneous broadening. As a result, an echo is discernible in in the time-domain data of Fig. 9b, despite the fact that is essentially completely masked in the individual contributions of Fig. 7 (right column).

During the fit, the absolute intensity of the experimental signal was not considered. Likewise, the proportionality constants in Eq. 4, $d\mu/dq$ and $d\alpha/dq$, are not known. Therefore, the three parameters $\sigma_{\mu}$, $\sigma_{\alpha}$, and $\sigma_{\mathrm{M}}$ are only known modulo an overall scaling factor, and only their relative contributions to the signal (e.g.., $\sigma_{\mu}/\sigma_{\mathrm{M}}$ and $\sigma_{\alpha}/\sigma_{\mathrm{M}}$) can be extracted. The sign of the overall signal can however be determined. The sign of the overall signal is given by the signs of $\sigma_{\mathrm{M}}$ times that of $d\alpha/dq$ (since $d\mu/dq$ enters quadratically in Eq. 10). To reproduce the experimentally observed sign, we either have $\sigma_{\mathrm{M}}$ negative and $d\alpha/dq$ positive (which is what we assumed here), or vice versa.

IV Discussion and Conclusion

We have shown that the experimental 2D-Raman-THz response of liquid water can be explained almost quantitatively by describing the hydrogen bond bending modes with the help of a very simple anharmonic oscillator model. The model contains only six independent parameters, i.e., $\sigma_{\mu}/\sigma_{\mathrm{M}}$ and $\sigma_{\alpha}/\sigma_{\mathrm{M}}$ as well as $\omega_{0}$, $\Delta\omega$, $\tau_{c}$, and $T_{1}$. Despite the fact that the instrument response function masks the information content of the molecular 2D response to a significant extent, it is still sufficient to determine these parameters with confidence (Fig. 10), if we assume the response originates from predominantly one of the two bands in the observation window of the experiment. We tentatively conclude that this dominating band is the hydrogen-bond stretching vibration at 170 cm*-1*, but we cannot exclude additional contributions from the hydrogen-bond bending vibration around 50 cm*-1*, or from couplings between both bands. The librational mode around 600 cm*-1*, on the other hand, is completely outside our experimental observation window (Fig. 5).

The physical interpretation of $\sigma_{\mu}$, $\sigma_{\alpha}$, and $\sigma_{\mathrm{M}}$ is not straight forward, since the low frequency modes of water are collective in nature and presumably delocalized to a certain extent. It is nonetheless worth noting that the biggest contribution to the signal originates from mechanical anharmonicity $\sigma_{\mathrm{M}}$ and the dipole nonlinearity $\sigma_{\mu}$, while the contribution from the nonlinearity in the polarizability is negligible with $\sigma_{\alpha}=0.1\sigma_{\mathrm{M}}$ (Tab. 2). This observation is in stark contrast to results from recent MD work, which revealed an echo for the TRT pulse sequence.Hamm2012a ; Ito2015 ; Ikeda2015 That echo originated from the hindered rotation band around 600 cm*-1*, which we however do not observe in the experiment due to the limited bandwidth. A large nonlinearity in the polarizability for that mode can be understood from the fact that a strict $\Delta J=2$ selection rule would apply for a Raman interaction, together with a $\Delta J=1$ selection rule for a THz interaction, in the limiting case of a free rather than hindered rotor.Hamm2017 This explanation obviously does not apply for the hydrogen-bond bending and stretching vibrations.

Regarding the line shape function, on the other hand, the parameters for $\tau_{c}$, $\Delta\omega$ and $T_{1}$ quantify the degree of inhomogeneous broadening of the hydrogen-bond-bending vibration (Tab. 2). That is, with $\Delta\omega\tau_{c}\approx 5$, the lineshape function is in the “slow-modulation”, or “quasi-inhomogeneous” limit.hamm11 At the same time, vibrational relaxation with $T_{1}=$250 fs contributes only 20 cm*-1* to the total linewidth. The correlation time $\tau_{c}=$370 fs, in turn, is a measure of the lifetime of the hydrogen-bond networks giving rise to the hydrogen-bond stretching vibration. The typical lifetime of a single hydrogen bond is 1 ps,asbury04a ; yere03 ; eaves05 ; perakis2011 , hence we conclude that those modes are delocalized over $\approx 3$ hydrogen bonds.

The quality of the fit of Fig. 9 is much better than that of much more sophisticated calculations based on a water force field in connection with MD simulations.Hamm2014 The problem with the MD approach arises from the fact that a water force field needs to describe the thermodynamics and dynamics of water reasonably well in the first place, while the anharmonicities, that are the bottleneck of the 2D-Raman-THz signal, result from these constraints only in an indirect way. Anharmonicities are typically not fitted explicitly, except if very specific effects, such as nuclear quantum effects,Habershon2009 are to be described. Furthermore, in particular when polarizability is included in a water force field, the number of parameters is typically very large, and the problem is underdetermined. Electrical anharmonicity has a lot to do with the redistribution of charges during the motion of water molecules, and being able to quantify it might reveal guidelines to design better water models.Sidler2018 To that end, it would be important to infer the anharmonicity parameters $\sigma_{\mu}$, $\sigma_{\alpha}$, and $\sigma_{\mathrm{M}}$ from a MD simulation of a realistic water force field.

In conclusion, we think that 2D-Raman-THz spectroscopy has the highest information content as to date with regard to the intermolecular degrees of freedom water, but learning how to extract that information from the experimental response is challenging. The present work constitutes a significant step in that direction.

Acknowledgments: We thank Yoshitaka Tanimura for very valuable discussions. The work has been supported by the Swiss National Science Foundation (SNF) through the NCCR MUST.

Appendix A Position Operator in an Anharmonic Eigenstate Basis

In the cubic anharmonic oscillator model, the harmonic oscillator Hamiltonian $\hat{H_{0}}$ is perturbed by $\sigma_{M}\hbar\omega q^{3}$.

[TABLE]

where $\hat{q}=\sqrt{m\omega/\hbar}\hat{x}$ is the unitless position operator, and $\sigma_{M}$ the size of the perturbation. To first order in $\sigma_{M}$, the eigenfunctions of $\hat{H}$ are expressed as linear combination of the harmonic eigenfunctions $|\varphi_{n}\rangle$. Claude1991

[TABLE]

with

[TABLE]

Mechanical anharmonicity does not affect the position operator $\hat{q}$ per se, but rather changes its matrix representation due to a change of the basis functions $\{|\varphi_{i}\rangle\}\rightarrow\{|\varphi^{(anh)}_{i}\rangle\}$. We will call $\hat{q}$ in a harmonic eigenfunction basis $\boldsymbol{q_{H}}$, and $\boldsymbol{q_{A}}$ if it is expressed in anharmonic eigenfunctions. Expressing $(\boldsymbol{q_{A}})_{i,j}\equiv\langle\varphi^{(anh)}_{i}|\hat{q}|\varphi^{(anh)}_{j}\rangle$ in terms of $(\boldsymbol{q_{H}})_{k,l}\equiv\langle\varphi_{k}|\hat{q}|\varphi_{l}\rangle$, using Eq. 23, we get a large collection of terms:

[TABLE]

where terms higher than first order in $\sigma_{M}$ have been discarded. The matrix form of the position operator in harmonic oscillator eigenbasis is (see Eq. 3):

[TABLE]

A matrix with $(\boldsymbol{q_{H}})_{i,j+1}$ is the matrix $\boldsymbol{q_{H}}$ shifted up by one row:

[TABLE]

Hence, the elements $\sqrt{j}\delta_{i+1,j}$ are transformed to $\sqrt{i+1}\delta_{i,j}$ terms in that shifted matrix, and $\sqrt{i}\delta_{i-1,j}$ of $\boldsymbol{q_{H}}$ to $\sqrt{i}\delta_{i-2,j}$. Using this procedure, the shifted matrices in Eq. 25 become:

[TABLE]

which include zero- ($\delta_{i,j}$), two- ($\delta_{i-2,j}$ and $\delta_{i,j-2}$), and four-quantum ($\delta_{i-4,j}$ and $\delta_{i,j-4}$) contributions. Collecting all terms, we obtain:

[TABLE]

We see that the $\delta_{i-4,j}$ elements cancel (likewise for $\delta_{i,j-4}$), and only zero- and two-quantum transitions remain in $\boldsymbol{q_{A}}$, in addition to the one quantum transitions from $\boldsymbol{q_{H}}$. Inserting Eq. LABEL:eq:fourQ into Eq. 25 we get

[TABLE]

where the second term describes the correction to the position matrix due to anharmonicity with the final expression of $\boldsymbol{q_{M}}$ given in Eq. 8.

Appendix B Temperature Independence of the Response Function

Temperature enters via the thermal population of the initial density matrix $\rho_{\mathrm{eq}}$. In IR spectroscopy, one usually assumes that only the ground state is initially populated due to $\hbar\omega\gg k_{B}T$, but that assumption no longer holds in THz spectroscopy, and the contributions from many initial states have to be considered. Fig. 12 exemplifies that for $R_{\mathrm{RTT}}^{\mu}$ by showing all contributing Feynman diagrams sorted according to their initial state. Since that response functions contain all four nonzero peaks (a), (b), (c), and (d), the result is universal and applies in the same way for the other response functions. It can be seen that the overall intensities do not depend on initial state and hence also not on temperature. However, a prerequisite for this to work is the assumption that the lineshape functions depend on energy level differences only (Eq. 14), and that the energy spacing between states is equidistant, i.e. that of a harmonic oscillator, and that the effect of mechanical anharmonicity only enters via the softened selection rules that allow for zero- and two-quantum transitions. Fig. 3 demonstrates that this is a good approximation (see discussion in the main text).