Interpretation of the THz-THz-Raman Spectrum of Bromoform

TL;DR

This study interprets the THz-THz-Raman spectrum of bromoform using an improved RDM model, simplifying the analysis by a self-consistent time reference and eliminating the need for large multi-quanta dipole elements.

Contribution

The paper introduces a self-consistent time coordinate framework and an analytical parametrization of the RDM model for better spectrum interpretation.

Findings

Simplified spectrum analysis with a self-consistent time reference.

Elimination of the need for large multi-quanta dipole matrix elements.

Agreement of RDM parameters with ab initio calculations.

Abstract

Nonlinear THz-THz-Raman (TTR) liquid spectroscopy offers new possibilities for studying and understanding condensed-phase chemical dynamics. Although TTR spectra carry rich information about the systems under study, the response is encoded in a three-point correlation function comprising of both dipole and polarizability elements. Theoretical methods are necessary for the interpretation of the experimental results. In this work, we study the liquid-phase dynamics of bromoform, a polarizable molecule with a strong TTR response. Previous work based on reduced density matrix (RDM) simulations suggests that unusually large multi-quanta dipole matrix elements are needed to understand the measured spectrum of bromoform. Here, we demonstrate that a self-consistent definition of the time coordinates with respect to the reference pulse leads to a simplified experimental spectrum. Furthermore, we analytically derive a parametrization for the RDM model by integrating the dipole and polarizability elements to the 4th order in the normal modes, and we enforce inversion symmetry in the calculations by numerically cancelling the components of the response that are even with respect to the field. The resulting analysis eliminates the need to invoke large multi-quanta dipole matrix elements to fit the experimental spectrum; instead, the experimental spectrum is recovered using RDM simulations with dipole matrix parameters that are in agreement with independent ab initio calculations. The fundamental interpretation of the TTR signatures in terms of coupled intramolecular vibrational modes remains unchanged from the previous work.

Figures

Tables

Table 1. Table 1: Values of the transition dipole matrix, obtained from ab initio (CCSD/aug-cc-pVTZ) calculations 31 and via fitting of the RDM models in Ref.  31 and in the current work.

$⟨i|\hat{M}|j⟩$	#Q	CCSD31	RDM (Ref. 31)	RDM (current)
$|000⟩\to |100⟩$	1	1.00	$1.00$ (fixed)	$1.00$ (fixed)
$|000⟩\to |001⟩$	1	1.01	$1.60$ (fixed)	$1.60$ (fixed)
$|100⟩\to |001⟩$	2	0.04	$0.009\left({\mu }_4\right)$	$0.000\left({\mu }_4\right)$
$|200⟩\to |101⟩$	2	0.06	$0.000\left({\mu }_4^{\prime }\right)$	$0.000\left(\sqrt{2}{\mu }_4\right)$
$|101⟩\to |002⟩$	2	0.06	$0.000\left({\mu }_4^{\prime \prime }\right)$	$0.000\left(\sqrt{2}{\mu }_4\right)$
$|100⟩\to |002⟩$	3	0.00	$0.364\left({\mu }_8\right)$	$-0.014\left(\frac{\sqrt{2}}{2}{\mu }_8\right)$
$|001⟩\to |200⟩$	3	0.01	$0.442\left({\mu }_7\right)$	$0.013\left(\frac{\sqrt{2}}{2}{\mu }_7\right)$
$|101⟩\to |300⟩$	3	0.02	$0.034\left({\mu }_7^{\prime }\right)$	$0.023\left(\frac{\sqrt{6}}{2}{\mu }_7\right)$
$|001⟩\to |300⟩$	4	-	$0.023\left({\mu }_{11}\right)$	$0.000\left(\frac{\sqrt{6}}{6}{\mu }_{11}\right)$
$|200⟩\to |002⟩$	4	-	$0.092\left({\mu }_{12}\right)$	$-0.005\left(\frac{1}{2}{\mu }_{12}\right)$

Table 2. Table 2: Hamiltonian and transition polarizability parameters of the RDM model.

$\hslash {\omega }_1$	4.7(fixed)	$\hslash {\omega }_3$	6.6(fixed)
${\Delta }_1/\hslash {\omega }_1$	-0.52%	${\Delta }_3/\hslash {\omega }_3$	-0.10%
${\beta }_{12}/\hslash {\omega }_1$	-0.05%	${\beta }_{13}/\hslash {\omega }_3$	+0.14%
${\alpha }_1$	1.00 (fixed)	${\alpha }_2$	-0.003

Equations

$$\frac{\partial \rho (t)}{\partial t}=-\frac{i}{\hslash }\left[H(t),\rho (t)\right]-\Gamma \rho (t),$$

$$H(t;t_1)=H-M\cdot E_T(t)-M\cdot E_T(t-t_1),$$

$$\rho (t+dt;t_1)$$

$$-2\Gamma \rho (t-dt;t_1)dt,$$

$$S^{Th}(t_2;t_1)=\Re \left(\mathrm{Tr}(\Pi \cdot \rho (t_2;t_1)P)\right),$$

$$P_{nm}=i,P_{mn}=-i,\text{and}{P}_{nn}=1.$$

$${\tilde{S}}^{Th}(f_2;f_1)=|\mathrm{FT}(S^{Th}(t_2;t_1)|$$

$$\delta =\iint d{f}^2{\left[log({\tilde{S}}^{Ex}(f_2;f_1))-log({\tilde{S}}^{Th}(f_2;f_1))\right]}^2$$

$$\Delta =\delta +\lambda (h_i^2+{\mu }_j^2+{\alpha }_k^2),$$

$$(t_2^{\prime },t_1^{\prime })=\{\begin{array}{ll}(t_2,t_1),& \text{region A}\\ (t_2+t_1,t_1),& \text{region B}\end{array}$$

$$(t_2^{\prime \prime },t_1^{\prime \prime })=\{\begin{array}{ll}(t_2^{\prime },t_1^{\prime })=(t_2,t_1),& \text{region A}\\ (t_2^{\prime },-t_1^{\prime })=(t_2+t_1,-t_1),& \text{region B}\end{array}$$

$$S(E)={\epsilon }_{1}E+{\epsilon }_2{E}^2+{\epsilon }_3{E}^3+{\epsilon }_4{E}^4+\dots ,$$

$$S(-E_T,-E_T,-E_{NIR})=-S(E_T,E_T,E_{NIR}),$$

$${\epsilon }_2={\epsilon }_4=\dots =0.$$

$$H_2(t;t_1)=H-M\cdot E_T(t)\text{and}$$

$$H_1(t;t_1)=H-M\cdot E_T(t-t_1).$$

$${\sigma }_{3^{rd}}\propto ⟨10|\frac{{\partial }^{2}M}{\partial Q_1\partial Q_2}|01⟩E_T.$$

$${\sigma }_{4^{th}}\propto ⟨10|\frac{\partial M}{\partial Q_2}|11⟩⟨11|\frac{\partial M}{\partial Q_1}|01⟩E_T^2.$$

$$\frac{{\partial }^{2}M}{\partial Q_1\partial Q_2}\lesssim E_T\frac{\partial M}{\partial Q_2}\frac{\partial M}{\partial Q_1}.$$

$$S^{Th}=S^{Th}(E_T,E_T)+S^{Th}(-E_T,-E_T).$$

$$\hat{H}$$

$$\hat{M}$$

$$+{\hat{M}}_{3q}({\mu }_7,{\mu }_7^{\prime },{\mu }_8)+{\hat{M}}_{4q}({\mu }_{11},{\mu }_{12}),$$

$${\hat{Q}}_i=\sqrt{\frac{2m_i{\omega }_i}{\hslash }}{\hat{q}}_i={\hat{a}}_i^{†}+{\hat{a}}_i.$$

$$\hat{H}$$

$$+{\Delta }_1\left({\hat{a}}_1^{†}{}^2{\hat{a}}_1^2+{\hat{a}}_2^{†}{}^2{\hat{a}}_2^2\right)+{\Delta }_3{\hat{a}}_3^{†}{}^2{\hat{a}}_3^2,$$

$$\hat{M}$$

$$+\frac{1}{3!}{\mu }_6\left({\hat{Q}}_1^3+{\hat{Q}}_2^3\right)+\frac{1}{2!}{\mu }_6^{\ast }{\hat{Q}}_1{\hat{Q}}_2\left({\hat{Q}}_1+{\hat{Q}}_2\right)+\frac{1}{2!}{\mu }_7\left({\hat{Q}}_1^2+{\hat{Q}}_2^2\right){\hat{Q}}_3+{\mu }_7^{\ast }{\hat{Q}}_1{\hat{Q}}_2{\hat{Q}}_3+\frac{1}{2!}{\mu }_8\left({\hat{Q}}_1+{\hat{Q}}_2\right){\hat{Q}}_3^2$$

$$+\frac{1}{3!}{\mu }_{10}^{\ast }{\hat{Q}}_1{\hat{Q}}_2\left({\hat{Q}}_1^2+{\hat{Q}}_2^2\right)+\frac{1}{2!2!}{\mu }_{10}^{\ast \ast }{\hat{Q}}_1^2{\hat{Q}}_2^2+\frac{1}{3!}{\mu }_{11}\left({\hat{Q}}_1^3+{\hat{Q}}_2^3\right){\hat{Q}}_3+\frac{1}{2!}{\mu }_{11}^{\ast }{\hat{Q}}_1{\hat{Q}}_2\left({\hat{Q}}_1+{\hat{Q}}_2\right){\hat{Q}}_3$$

$$+\frac{1}{2!2!}{\mu }_{12}\left({\hat{Q}}_1^2+{\hat{Q}}_2^2\right){\hat{Q}}_3^2+\frac{1}{2!}{\mu }_{12}^{\ast }{\hat{Q}}_1{\hat{Q}}_2{\hat{Q}}_3^2.$$

$${\hat{Q}}_i^n$$

$$+\dots +\sum _{perm}{\hat{a}}_i^{†}{}^2{\hat{a}}_i^{(n-2)}+\sum _{perm}{\hat{a}}_i^{†}{\hat{a}}_i^{(n-1)}+{\hat{a}}_i^n$$

$$H^{\prime }=D^{†}HD$$

$$M^{\prime }=D^{†}MD,{\Pi }^{\prime }=D^{†}\Pi D.$$

Full text

1 Introduction

Two dimensional time resolved spectroscopy encompasses techniques in which a sequence of three or four coherent laser pulses is used to measure the response of the system with respect to the time delays. The direct Fourier transform of this response corresponds to a frequency-frequency correlation function that carries rich information about the anharmonicity, coupling and non-linearity of the vibrational and electronic states. These characteristics elucidate the underlying dynamics and various energy dissipation processes, such as homogeneous and inhomogeneous broadening.

Each laser pulse in the sequence creates or destroys coherences via one of two physical processes: optical absorption (1st order in the field interaction) or Raman scattering (2nd order in the field interaction). These processes can take place in the visible (VIS) regime probing electronic states, the infrared (IR) regime measuring high frequency intra-molecular vibrational modes, or the terahertz (THz) regime probing low-frequency modes such as lattice phonons or inter-molecular vibrations. Combination of these various laser pulses leads to a range of different spectroscopic methods such as 2D-IR, 2D-VIS, VIS-IR, 2D-Raman, etc., each with its own applicability.

2D-IR spectroscopy was one of the first methods used to investigate the vibrational properties of liquids, and it has helped address topics that include structural fluctuations in water 1, 2, 3, 4, proton shuttling in acidic solutions 5, 6 and protein dynamics 7, 8, 9. While 2D-IR has proven useful for studying high-frequency vibrational modes, there is strong motivation to extend multidimensional spectroscopy to the low-frequency regime. Important physical chemical processes, such as collective solvent motions 10, 11, 12 and dynamics of large biomolecules13, 14, 15 are driven by processes that takes place in the THz range. The first method developed to study liquids in the THz domain was 2D-Raman spectroscopy16. This method is 5th order in the laser field, and it has been shown to suffer from cascading effects where the 5th order response is plagued by contributions from higher intensity 3rd order processes17. It is only with difficulty that these cascaded processes were overcome to yield the true 5th order 2D Raman response.18, 19

The advent of powerful THz sources 20, 21 has enabled 3rd order hybrid spectroscopic techniques, in which one or two of the 2nd order Raman processes are replaced by 1st order THz absorption: RTT, TRT, and TTR 22. The RTT and TRT methods measure a THz emission and were first employed to study water and aqueous solutions 23, 12, 24, 25. TTR builds upon the pulse-detection methods developed in THz-Kerr-effect spectroscopy,26, 27 which has been successfully applied in both polar and non-polar liquids 28, 29. The TTR approach was recently developed in our group 30, 31 and used to measure the 3rd order nonlinear response of bromoform 32 (CHBr3), carbon tetrachloride (CCl4) and dibromodichloromethane (CBr2Cl2). These halogen liquids are ideal test systems for TTR because they exhibit heavy intra-molecular modes in the THz frequency regime and are Raman-active due to their strong polarizablity.33

All three Raman-THz hybrid methods (i.e., RTT, TRT, and TTR) are complementary and measure different correlation functions involving the dipole and polarizability surfaces. These responses carry rich information about the systems under study, including the nonlinearity of the dipole and polarizability surfaces and the anharmonicity and mechanical coupling of the various vibrational motions. However, this information is encoded in complex, three-point correlation functions that must be disentangled with the aid of theoretical and computational methods 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 30, 31, 48. This problem of interpretation is further complicated by the number of vibrational states that need to be considered. Whereas infrared-active modes are usually in their vibrational ground state at room temperature and typically only involve single excitations upon illumination, THz-active modes can be thermally excited at room temperature, and multiple transitions between the different states are possible.

In our previous work 31, we developed a reduced density matrix (RDM) model to understand the TTR spectrum in liquid bromoform. When the parameters for the nonlinear dipole elements in the RDM model were fit to best reproduce the experimental TTR signal, it was found that they assumed unexpectedly large values that did not agree with our accompanying ab initio electronic structure calculations. Here, we reconcile this apparent inconsistency by developing a more complete description of the TTR spectra from both the theoretical and experimental perspective.

2 Method Development

We begin by reviewing the RDM model of Ref. 30, 31. In section Interpreting the Experimental Data (2.1), we describe the time-coordinate transformation necessary to reinterpret the experimental spectrum. In sections 2.2 and 2.3, we explain the development of a new RDM model that fully accounts for the liquid symmetry and the symmetries of the dipole and polarizability matrices to 4th order in the normal modes.

In RDM, we propagate the reduced density matrix $\rho(t)$ for a single bromoform molecule by solving the Liouville–von Neumann equation 49,

[TABLE]

where $\Gamma$ is a constant phenomenological relaxation matrix that accounts for the interaction with the bath by population relaxation ($\Gamma_{nn}=\frac{1}{\tau_{1}}$, where $\tau_{1}=100$ fs) and coherence dephasing ($\Gamma_{nm}=\frac{1}{\tau_{2}}$, where $\tau_{2}=2$ ps). We tested longer population relaxation times up to 10 ns, and we found no change in the computed TTR response.

The time-varying Hamiltonian, $H(t)$, describes interactions of the system with the THz fields,

[TABLE]

where $E_{T}(t)$ is the experimentally measured pulse shape 31, $H$ is the system Hamiltonian, and $M$ is the transition dipole matrix. Using $E_{T}(t)$ in the simulation, instead of a simple delta function $E\delta(t)$ ensures that the computed spectrum is in fact the molecular response convoluted with the instrument response function (IRF) and can be directly compared to the raw experimental data.

We propagate equation 1 forward in time by numerical integration, noting that the following mixed central-forward scheme provides good numerical stability:

[TABLE]

where the commutator is discretized by central difference, which preserves time reversibility, and the phenomenological dissipation is discretized by forward difference. We then calculate the TTR signal as the field emitted by the final Raman process,

[TABLE]

where $\Pi$ is the transition polarizability, another important model parameter, and $P$ is a matrix with elements

[TABLE]

$P$ ensures that the signal is emitted $90^{\circ}$ out-of-phase with respect to the induced polarization 49.

We compute the frequency response from the absolute value of the Fourier transform, normalized to its maximum value,

[TABLE]

For a given realization of parameters in $H$, $M$, and $P$, comparison of the simulated and experimental TTR signals is made via the logarithmic error function:

[TABLE]

We fit the RDM model by minimizing $\delta$ with respect to $H$, $M$, and $P$, subject to $L_{2}$-norm regularization. Specifically, we minimize

[TABLE]

where $h_{i}$, $\mu_{i}$ and $\alpha_{i}$ are the parameters of the Hamiltonian $H=H(h_{i})$, transition dipole $M=M(\mu_{i})$ and transition polarizability $\Pi=\Pi(\alpha_{i})$ operators.

We note that Eq. 1 only accounts for Markovian dissipative processes. This description may of course be extended to include more complex relaxation phenomena such as coherence and population transfer in the context of non-Markovian dynamics,50, 51, 52 which have been shown to be important in some liquids.

2.1 Interpreting the Experimental Data

The experimental setup30 consists of two THz pump pulses, THz1 and THz2, which are polarized along the Y and X directions and separated by a time delay, $t_{1}=t_{\text{THz}_{1}}-t_{\text{THz}_{2}}$. A near-infrared (NIR) optical probe, polarized along X, induces a Y-polarized Raman response from the system at time delay $t_{2}=t_{\text{NIR}}-t_{\text{THz}_{1}}$. Due to the symmetry of the third-order response, a Raman signal is only detected after the second THz pump has interacted with the system. Experimentally, $t_{1}$ is scanned by changing the path length that the THz2 field travels, while the THz1 field remains fixed in time. A consequence of this design is that the third-order response will occur at $t_{2}=0$ when THz1 is the second field interaction ($t_{1}>0$), but at a time $t_{2}=t_{1}$ when THz2 is the second interaction ($t_{1}<0$). Qualitatively, a diagonally skewed TTR response is observed in the $t_{1}<0$ region. This must be properly accounted for in the interpretation of the spectrum.

To illustrate the importance of correctly including this skew in the analysis, we consider how a Liouville pathway in a generic system composed of two vibrational modes, changes when evaluated at $t_{1}>0$ and $t_{1}<0$. In Fig. 1(a), we label these two regions A and B, corresponding to pulse ordering THz2-THz1-NIR ($t_{1}>0$) and THz1-THz2-NIR ($t_{1}<0$), respectively. When the pathway is sampled in region A of the response, the first THz pulse excites coherence $\lvert 10\rangle\langle 00\rvert$, which oscillates at frequency $\omega_{1}$ and acquires phase $e^{i\omega_{1}t_{1}}$, while the second THz pulse switches the coherence to $\lvert 01\rangle\langle 00\rvert$ which oscillates at frequency $\omega_{2}$ and acquires phase $e^{i\omega_{2}t_{2}}$. Overall, this physical process generates a signal that is proportional to $e^{i\omega_{1}t_{1}}e^{i\omega_{2}t_{2}}$. The same pathway, when sampled in region B, generates a signal that is proportional to $e^{-i\omega_{1}t_{1}}e^{i\omega_{2}(t_{2}+t_{1})}$, or rearranging, $e^{i(\omega_{2}-\omega_{1})t_{1}}e^{i\omega_{2}t_{2}}$. Hence, taking the Fourier transform of the full $t_{1}$ domain (including both regions A and B) results in two distinct peaks at $(\omega_{2},\omega_{1})$ and $(\omega_{2},\omega_{2}-\omega_{1})=(-\omega_{2},\omega_{1}-\omega_{2})$ that otherwise represent the same physical process. This is undesirable and creates the opportunity for misinterpreting the TTR spectrum.

To ensure that each physical process is represented by a single distinct peak, we perform a coordinate transformation of the time response before computing the frequency response by FT. This transformation consists of skewing the response along $t_{2}$,

[TABLE]

and then flipping it with respect to $t_{1}$,

[TABLE]

as shown in Fig. 1(c). With these transformed time coordinates, our illustrative Liouville pathway acquires the same phase factors $e^{i\omega_{1}t_{1}^{\prime\prime}}e^{i\omega_{2}t_{2}^{\prime\prime}}$ when sampled in region B as in region A. This transformation simplifies the resulting TTR spectrum and eliminates the appearance of redundant peaks.

2.2 Enforcing Inversion Symmetry in the RDM Model

Liquids are isotopic, such that any response function,

[TABLE]

must obey inversion symmetry with respect to an applied electric field, $E$. The TTR response must also obey this symmetry,

[TABLE]

and requires that even-order contributions to the experimentally observed response must vanish,

[TABLE]

This means that only scattering processes that involve an odd number of photons, not counting the final emission, can generate an experimental TTR signal. In the case of liquid bromoform, the final Raman process is linear and it only involves a one-photon absorption 31. Therefore, the total number of THz photon interactions must be even and the leading contribution to the TTR signal is 3rd in the field, involving the absorption of two THz photons and one NIR photon. In other words, the TTR response scales linearly with each of the THz fields.

Fig. 2(b) lists all possible 3rd order interactions in which two THz photons are involved: 1-1-1, 2-0-1 and 0-2-1. The TTR experiment further isolates the desired 1-1-1 response by employing differential chopping 30. In the RDM simulations, we implement differential chopping by separating the 1D responses $S_{2}^{Th}$ and $S_{1}^{Th}$ from the total response $S^{Th}$, as illustrated in Fig. 2(a). We compute these individual responses by propagating the dynamics under the partial Hamiltonians

[TABLE]

This simulation procedure removes the unwanted field interactions, shaded purple in Fig. 2, from the 3rd response. However, the RDM simulated response 30, 31 also includes 4th order interactions, unlike the experiment where even order processes vanish due to the symmetry of the liquid, as explained above.

It is found that for the simulated TTR response, 4th order contributions can be comparable in magnitude to 3rd order ones, because the latter require at least one nonlinear dipole interaction222To prove that all 3rd order TTR pathways require a nonlinear dipole interaction, consider a 3rd order Liouville pathway which starts from a generic population state $\lvert a\rangle\langle a\rvert$. Interaction of the system with the first THz pulse changes the excitation state of the bra or ket of the density matrix by $n$ quanta to a $\lvert a+n\rangle\langle a\rvert$ coherence. Then, upon a second THz interaction, the system changes to either $\lvert a+n+m\rangle\langle a\rvert$ or $\lvert a+n\rangle\langle a+m\rvert$. Since the final Raman interaction only changes the system by a single quantum 31, and since all pathways must end in a population state, it follows that either $\lvert n+m\rvert=1$ or $\lvert n-m\rvert=1$. All paths involving $n=0$ or $m=0$ are removed by differential chopping, and therefore all surviving diagrams have both nonzero $n$ and $m$. As a result, at least one of the THz processes exchanges more than one quantum with the systems, and it is strictly nonlinear. while the former do not. For instance, in Fig. 2(e), we compare two similar pathways of 3rd and 4th order. In the 3rd order process, the two-quanta excitation from $\lvert 10\rangle$ to $\lvert 01\rangle$ is achieved by a single photon interacting with a nonlinear dipole element and has a scattering amplitude proportional to

[TABLE]

Meanwhile, in the 4th order process, the same excitation is achieved by two consecutive photons acting on the linear dipole via an intermediate state,

[TABLE]

For large THz fields, this “all-linear” 4th order can become larger than the 3rd order if the system has small dipole nonlinearities, i.e.,

[TABLE]

To address this issue, 4th order contributions to the simulated RDM response must be explicitly removed. We achieve this by simulating both the response to positive THz fields and negative THz fields, and then summing the two contributions as illustrated in Fig. 2(d),

[TABLE]

This procedure ensures that the final response function is odd with respect to the overall field and that the 4th order contributions, shaded green in Fig. 2(c), vanish. As a result, the simulated TTR spectrum correctly reflects the inversion symmetry of the liquid state and can be compared directly to the experimental spectrum.

2.3 The Hamiltonian, Dipole, and Polarizability Surfaces

The bromoform molecule has three vibrational degrees of freedom in the terahertz regime: two degenerate C-Br bending modes at 4.7 THz, $Q_{1}$ and $Q_{2}$, and a symmetric umbrella mode at 6.6 THz, $Q_{3}$ 53. Fig. 3(a) illustrates these motions and provides an energy-level diagram with the number of quanta in each mode indicated, using notation $\lvert Q_{1}Q_{2}Q_{3}\rangle$. In typical 2D-IR spectroscopy applications at room temperature, only the ground and first-excited state of each mode are involved. TTR spectroscopy probes modes that are low in energy, such that even at room temperature ($kT\approx 6.2$ THz), it is necessary to consider thermally accessible excited states. In the current study, the bandwidth of the THz pulse covers $BW\approx 8$ THz of energy, such that all the states below $kT+BW\lesssim 14.2$ THz should be considered. The highest energy state that is included in our model is the quadruplet $Q_{1}\otimes Q_{2}$, amounting to three quanta of energy in the degenerate mode ($3\times 4.7$ THz = $14.1$ THz).

The calculated TTR spectrum depends on the parameterization of $H=H(h_{i})$, $M=M(\mu_{i})$ and $\Pi=\Pi(\alpha_{i})$. In the RDM model of Ref. 31, we employed a harmonic Hamiltonian,

[TABLE]

the parameters for which were fixed to the linear absorption experimental spectrum, and a transition dipole of the form

[TABLE]

where $\hat{M}_{nq}$ represent the nonlinear blocks of the dipole operator comprised of $n$-quanta $(n>1)$ transitions, which were determined by fitting to the nonlinear TTR experimental spectrum. The $\hat{Q}_{i}$ are non-dimensional normal modes which can be expressed in terms of the creation and annihilation operators $\hat{a}_{i}^{\dagger}$ and $\hat{a}_{i}$,

[TABLE]

Combining Eqs. 21 and 20 and integrating the linear part of the dipole in the $\lvert Q_{1}Q_{2}Q_{2}\rangle$ basis leads to the matrix parameterization of Ref. 31 (Fig. 3(b)).

In this work, we develop a more comprehensive model for the Hamiltonian, which includes the anharmonicity $\Delta_{i}$ of the single molecule normal modes, as well as mode coupling $\beta_{ij}$ originating from the condensed-phase environment. We employ a quantum-conserving Hamiltonian of the form {strip}

[TABLE]

and integrate it in the $\lvert Q_{1}Q_{2}Q_{2}\rangle$ basis to obtain the matrix representation shown in Fig. 3(c). Note that the couplings between the doubly degenerate modes $Q_{1,2}$ and mode $Q_{3}$ are equivalent and labeled $\beta_{13}$ in Eq. 22.

Furthermore, we Taylor expand the dipole operator to 4th order in the normal modes and enforce the $Q_{1,2}$ symmetry by grouping together terms which are invariant to label exchange (1$\leftrightarrow$2): {strip}

[TABLE]

Whereas the dipole form in Ref. 31 (Eq. 20) uses a normal mode expansion of the linear part and an empirical parameterization of the nonlinear part, the current work expands the full dipole operator in the normal modes to 4th order. As a result, the parameters of the current model are associated with the partial derivatives of $\hat{M}$ with respect to the $\hat{Q}_{i}$. We express the $n^{th}$ power of $\hat{Q}_{i}$ as creation-annihilation $n$-tuples,

[TABLE]

and we keep only the pure creation (first term) and pure annihilation (last term); all other terms account for transitions involving less than $n$ quanta. We then integrate Eq. 23 in the $\lvert Q_{1}Q_{2}Q_{2}\rangle$ basis and obtain the model illustrated in Fig. 3(d). The polarizability operator is parameterized in a similar way; however, it is only expanded to second power in $Q_{i}$, since the final Raman interaction is expected to be linear.31

To fit the simulated RDM spectrum to the experimental TTR result, we keep $\hbar\omega_{1}$, $\hbar\omega_{3}$, $\mu_{1}$, $\mu_{3}$, and $\alpha_{1}$ fixed, while varying the remaining parameters. We perform the RDM dynamics in the diagonal representation of the Hamiltonian. As such, we diagonalize $H$,

[TABLE]

and apply the same transformation to the dipole and polarizability,

[TABLE]

performing the RDM simulations with $H^{\prime},M^{\prime}$ and $P^{\prime}$.

3 Results and Discussions

The results are presented in two parts. We first present the revised experimental spectrum obtained through the transformation of the time coordinates.We then discuss the results of the RDM model introduced in this work.

3.1 Revised Experimental Spectrum

The raw experimental TTR response is comprised of an orientational molecular response with superimposed oscillations arising from vibrational coherences. As in previous work27, 30, we focus on the vibrational component of the response, so we use an exponential fit along $t_{2}$ to de-trend the orientational response from the data.

In Figs. 4(a)-(c), we demonstrate how the vibrational signal changes upon application of the time-coordinate transformation from section Interpreting the Experimental Data (2.1). The first step of the transformation involves a skew along $t_{2}$, which makes the response function appear symmetric with respect to the $t_{1}=0$ line. The transformed signal is symmetric because swapping the order of the THz pulses does not change the response type; in both cases, we measure a TTR signal. This is unlike a TRT experiment where inverting the pulse order results in a RTT sequence which measures a different correlation function altogether 23. The second time transformation involves a flip with respect to $t_{1}$, which due to the symmetry of the response about $t_{1}=0$ simply averages the A and B regions of the data and improves the signal to noise ratio.

Panels (d) and (e) of Fig. 4 show the absolute value Fourier transforms of the original and transformed time responses, respectively. The revised spectrum is significantly simplified (panel (e)). Whereas the original spectrum features six intense peaks, the revised spectrum shows three features that are familiar in the context of conventional two dimensional spectroscopy. The most intense peak forms at $f_{1}=\omega_{3}-\omega_{1}$ and reports on the coupling between modes $Q_{1}$ and $Q_{3}$ via a nonlinear dipole interaction. The equivalent mode in the negative $f_{2}<0$ quadrant is closely related and shows a clear rephasing pathway. Note that in the revised spectrum, the $f_{2}<0$ quadrant corresponds to true rephasing, whereas in the original spectrum the peaks involve a mix of rephasing and higher overtone non-rephasing pathways. Interestingly, the most intense true-rephasing feature appears tilted at an angle; this tilt could be related to the degree of inhomogeneous broadening. Finally, in the non-rephasing quadrant along the first diagonal, the higher frequency features appear symmetric to those at lower frequency.

3.2 Parameterization of the RDM Model

We fit the parameters of the RDM model with respect to the revised experimental spectrum using the Levenberg-Marquardt algorithm implemented in Octave 54, 55, 56, minimizing the penalty function in equation 8. For each value of the regularization parameter $\lambda$, we start with random initial guesses for the RDM parameters, sampled uniformly from the $\left[-1,1\right]$ interval. The fitting problem is highly under-determined, featuring multiple local minima, due to the large number (order $10^{4}$) of Liouville paths that compete in the creation of a relatively small number of peaks. Additionally, the dipole and polarizability parameters can take negative values, which lead to peak cancellations. To address these challenges, we perform multiple minimizations, starting each time from a different random initial guess.

In Fig. 5(a) we plot all the converged solutions as function of the maximum fitted parameter $max(h_{i},\mu_{j},\alpha_{k})$, and color them according to the $L_{2}$ regularization strength, $\lambda$. Clustered to the left are under-fitted results with small fitting parameters but large errors, while to the right are the over-fitted results that match the spectrum at the expense of unphysical parameters. The optimal solution lies in the intermediate regime with $\lambda=100$, as indicated in Fig. 5(a). Three independent minimizations starting from different initial guesses converged to the same error with very similar parameters, suggesting relative robustness of this optimal paremeterization. We summarize the best-fit parameters in tables 1 and 2, which report the values of the transition dipole elements, and we show the full transition dipole matrix in Fig. 5(b).

The resulting RDM simulated spectrum shown in Fig. 5(c) is in excellent agreement with the revised experimental spectrum (Fig. 4(e)). This can be clearly seen in Fig. 5(d), which provides a side-by-side comparison between simulation and experiment along the three most important 1D traces from the 2D response.

It is found that use of the experimental pulse shape $E_{T}(t)$ in the RDM simulations is crucial for reproducing the experimental spectrum. Using simple $\delta$-functions to simulate the response, otherwise with the same RDM parameters, gives a radically different result, with just one dominant peak on the diagonal at $f_{2}=4.7$ THz.

The dipole nonlinearities fitted with the RDM model in this work are significantly smaller than those found in Ref. 31, and they agree well with results from ab initio electronic structure calculations (see table 1). Critically, the third-order fitted non-linearities in the current RDM model have the same order of magnitude as those from the ab initio electronic structure calculations, in contrast with the third-order terms from the fitted RDM model of Ref. 31. This central result demonstrates that the refined treatment of the RDM model and the experimental spectrum in the current work removes the central inconsistency of Ref. 31. Likewise, the fourth-order terms in the fitted RMD model of the current work are vastly smaller than those that were necessarily invoke in the previous study31. For the second-order terms, both the current and earlier RDM models have very small fitted values. It is found that the difference between the fitted and ab initio values for the second-order terms is insignificant. We confirmed this by running the current RDM model with the ab initio values for $\mu_{4}$ and all other fitted parameters unchanged, and finding that the resulting simulated spectrum is unaffected.

Finally, in Table 2, we report the fit RDM parameters for the anharmonicity $\Delta_{i}$ of the normal modes, their mechanical coupling $\beta_{ij}$, and the first order polarizability elements $\alpha_{1}$ and $\alpha_{3}$. We find small anharmonicities and mechanical coupling elements, which points to relatively harmonic, weakly coupled normal modes in bromoform. This conclusion is in agreement with vibrational second-order perturbation theory (VPT2) calculations performed previously.31 Finally, the polarizability matrix elements obtained in this current work are in good agreement with those obtained in our previous work 31, with very small nonlinear elements indicating that the final Raman interaction is indeed a linear process.

It is worth noting that the simulated spectrum has low-intensity peaks at $f_{2}\simeq 6.6$ THz that bear resemblance to those in the experimental spectrum. Given weakness of these peaks, they make a small contribution to the error function in the fitting and are thus poorly resolved. Additionally, in some fits, we found that the features at $f_{2}\simeq 5.5$ THz might be connected to the mechanical coupling between $Q_{1,2}$ and $Q_{3}$. Future refinement of the RDM model and higher resolution experiments are needed to interpret these subtle features.

4 Conclusions

Our original TTR study of bromoform 31 revealed large dipole nonlinearities which were inconsistent with ab initio electronic structure calculations. Here, we refine both the experimental and theoretical description of this system, resolving the inconsistency. The key refinements are described below.

First, we revised the experimental spectrum by introducing a simplifying time-coordinate transformation. Due to the symmetry of the experimental setup, the time response must be skewed and flipped when the ordering of the two THz pulses is changed before taking the Fourier transform for a correct interpretation. The new spectrum reveals fewer features and true rephasing signals which are symmetric with their non-rephasing counterparts, as expected.

Second, we have developed an RDM model that analytically includes all dipole and polarizability matrix elements up to 4th order. We also include mechanical anharmonicity and mode-coupling in the model Hamiltonian. This refined RDM model preserves the symmetry of the bromoform molecule.

Lastly, we rigorously account for the inversion symmetry of the liquid by excluding 4th order field interactions from the response. We achieve this by simulating the response of the system to both positive and negative THz fields and combining the two results. The final response is cubic in the field, in agreement with the experimentally measured TTR response.

The revised experimental spectrum is used to fit the parameters of the updated RDM model, leading to good agreement. The fitted nonlinearities are orders of magnitude smaller than found in our previous work and agree well with ab initio electronic structure calculations. Regardless, the original conclusion that nonlinearities in the dipole surface of the intramolecular vibrations drive the TTR response of bromoform remains unchanged.

5 Acknowledgements

The authors thank Ralph Welsch, Ian Finneran, Matthew Welborn, and Philip Shushkov for helpful discussions, as well as Peter Hamm for a sharing a copy of their forthcoming manuscript. This work was supported by the National Science Foundation (Grant CHE-1665467) and the Office of Naval Research (Grant N00014-10-1-0884).