State-resolved infrared spectrum of the protonated water dimer: revisiting the characteristic proton transfer doublet peak

Abstract

The infrared (IR) spectra of protonated water clusters encode precise information on the dynamics and structure of the hydrated proton. However, the strong anharmonic coupling and quantum effects of these elusive species remain puzzling up to the present day. Here, we report unequivocal evidence that the interplay between the proton transfer and the water wagging motions in the protonated water dimer (Zundel ion) giving rise to the characteristic doublet peak is both more complex and more sensitive to subtle energetic changes than previously thought. In particular, hitherto overlooked low-intensity satellite peaks in the experimental spectrum are now unveiled and mechanistically assigned. Our findings rely on the comparison of IR spectra obtained using two highly accurate potential energy surfaces in conjunction with highly accurate state-resolved quantum simulations. We demonstrate that these high-accuracy simulations are important for providing definite assignments of the complex IR signals of fluxional molecules.

---

1 Introduction

Protonated water clusters play a key role in chemistry and biology, yet their quantum nature continues to puzzle scientists.^1–5 One of the smallest clusters, the protonated water dimer or Zundel ion, H(H₂O)₂⁺, is particularly important and interesting as it can be understood as one of the building blocks of acidic water,^6,7 where its infrared (IR) spectral signatures appear.^8–10 Its structure of a proton being “sandwiched” by two water units leads to an unconventional bonding,^7,11 serving as a limiting case of a characteristic pattern^6,12 that is observed in a variety of protonated clusters.⁷ Its IR spectrum is of particular importance, as it reveals information about chemical bonding and the quantum nature of molecular vibrations. However, the IR spectrum of the Zundel ion is complicated – as such and due to tagging effects.^3,13–16 Unexpectedly, its complexity even increases upon deuteration.^17,18 Most dominant in the spectrum is the famous doublet near 1000 cm⁻¹, whose vibrational character had not been understood initially.¹³ Quantum dynamics simulations of the Zundel ion allow us to understand the spectrum in full detail. However, they are difficult to perform for the following three reasons: (i) its dimensionality with 15 vibrational degrees of freedom, (ii) the non-obvious quantum entanglement of the vibrational zero-order states, and (iii) the fluxional nature of the protons in the Zundel ion with several equivalent minima, connected by only shallow potential barriers. When combined as in case of the Zundel ion, these difficulties severely challenge quantum dynamics simulations even today.

To date, the multiconfiguration time-dependent Hartree (MCTDH) method has been the only one that could fully reproduce the experimental IR spectrum (save for tagging effects) and it revealed for the first time the characteristics of the doublet.^14,15,19 According to MCTDH analysis, it consists of a Fermi resonance, where low-frequency, low-intensity wagging (pyramidalization) motions of the water units entangle with the proton transfer motion. The analysis has later been confirmed by other simulations.^16,20,21

While MCTDH wavefunction propagation enables the understanding of the Zundel IR spectrum, the analysis of the vibrational transitions is cumbersome because information on the vibrational structure is only indirectly accessible through observables and low-dimensional probability densities. Further, the spectral resolution of wavefunction propagations is limited by the propagation time and the longer the propagation time, the more complicated the propagated wavefunction becomes. It is desirable to being able to obtain the actual vibrational excited wavefunction that is responsible for a particular peak in the IR spectrum, but this is difficult due to the sheer number of excited eigenstates that are present in such a system of high anharmonicity and dimensionality.

Here, we present new simulations that not only give more accurate vibrational spectra via wavefunction propagations but also reveal the full eigenstate spectrum with up to 900 eigenstates to high accuracy for a vibrational energy up to ∼1900 cm⁻¹. Two highly accurate potential energy surfaces (PESs) and dipole moment surfaces (DMSs) are used, namely the pioneering permutationally-invariant-polynomial-based HBB surfaces²² and the more recent BBSM neural-network-based surfaces,^23,24 which so far has not been used for fully, rigorous quantum simulations. In the following, by referring to one of the used PESs we imply the usage of the corresponding DMS as well.

Our full-dimensional quantum dynamics simulations consist of recently introduced vibrational tree tensor network states (TTNS) methodology,²⁵ as well as new, more accurate and more compact representations of the PES.^26,27 The TTNS method is based on the density matrix renormalization group (DMRG),^28,29 shares the same wavefunction ansatz as the multilayer (ML) MCTDH method, enables both time-independent and time-dependent simulations^25,30,31 (denoted as ti-TTNS and td-TTNS), and, most importantly here, the time-independent simulations are computationally much less demanding than in ML-MCTDH.²⁵ This allows us to fully identify the nature of each peak in the spectrum by systematically characterising the eigenstates.

The key finding of this study is that the interplay between the wagging and the proton motions is more complex than previously thought in a critical way. In particular, we disclose that astonishingly subtle changes in the energetics of the Zundel ion can lead to very different features in the IR spectrum. This not only has important implications to PES development and the quantum dynamics of fluxional molecules in general, but it will also open the door toward understanding the chemistry of (micro-)solvated species.

Previous MCTDH simulations on the HBB PES revealed that the characteristic doublet contains contributions of two states, one with one quantum in the proton transfer motion, and another one with one quantum in the water–water stretch motion and two quanta in the wagging motions. Here, in contrast, our simulations on the same PES reveal that there are not two but three dominant eigenstates that are responsible for the “doublet” in the IR spectrum. Besides the known contributions from a state with two quanta in the wagging motions, there is an additional dominant contribution from a state having four quanta – only discovered with our novel TTNS-based quantum dynamics method – giving overall three dominant peaks in the IR spectrum. In contrast, and in accord with the established knowledge from the MCTDH simulations, the newer BBSM PES displays again only two peaks for the doublet whereas the state with four quanta in the wagging motion is blueshifted and can indeed be identified as satellite peak of the doublet in the experimental spectrum, in support of our assignment.

2 Setup

Our benchmark-quality simulations are based on the Hamiltonian in polyspherical coordinates set up by Vendrell et al.^15,19,32 We fitted both the HBB and the BBSM PESs into a form more suitable for (ML-)MCTDH and TTNS simulations. We use the recently developed Monte Carlo canonical polyadic decomposition (MCCPD) approach for refitting the PESs.^26,27 The refitting introduced a mean/root-mean-square error of ∼0.2 cm⁻¹/∼11 cm⁻¹. For comparison, the mean/root-mean-square error of the original cluster expansion of the HBB PES used in ref. 14, 15, 19 and 32 is about ∼24 cm⁻¹/∼4 × 10³ cm⁻¹.²⁷ To validate the accuracy of both the two PESs and, as additional validation, of the two MCCPDs of the respective PESs, we compare the energies of the PESs from samples of ∼2000 configurations with direct, “gold standard” basis-set-extrapolated coupled cluster singles and doubles with perturbative triples (CCSD(T)) electronic structure energies, see ESI† for details. With these ∼2000 configurations, for the HBB PES we obtain a mean absolute error (MAE) of 203 cm⁻¹, whereas for the BBSM PES we obtain an improved accuracy with an MAE of 75 cm⁻¹. The MCCPDs have essentially the same accuracy.

We represent the wavefunction as TTNS, which has the same mathematical structure as the wavefunction used in the ML-MCTDH method.^25,33 It differs, however, by the use of more efficient algorithms that are based on a generalization of the density matrix renormalization group (DMRG).^25,28,29 For the HBB /BBSM PES we computed up to ∼900/550 of the lowest eigenstates, despite the large density of states that is evident in this fluxional cluster. We verify the eigenstate “stick” spectrum by comparing with the spectrum obtained from wavefunction propagation using the time-dependent DMRG for TTNSs.^30,31 The Fourier-transformed propagation-based spectrum allows us to obtain the full spectrum with broadened peaks, while the stick spectrum reveals every single detail. Details of the employed methodology, the numerical parameters, convergence tests, as well as additional results are presented in the ESI.†

3 Results and discussion

Fig. 1 displays the experimental spectrum^13,17 and the spectra simulated using either of the two used PESs. For the HBB /BBSM PES we show in Fig. 1 up to ∼900/550 states up to ∼1900 cm⁻¹/∼1700 cm⁻¹; note that most of these states are not IR-active so that much fewer lines are visible in the spectra. The total spectrum, including IR inactive states, is shown as small lines on the abscissa. The higher the energy, the higher the density of states, thus indicating the large spectral complexity of the Zundel ion. The energies of the IR inactive states differ substantially between the two PESs, which will be discussed in more detail elsewhere. Compared to the propagation-based spectrum (td-TTNS), the stick spectrum (ti-TTNS) has slightly different intensities and the transitions are centered at slightly lower wavenumbers. This is because the stick spectrum is not convoluted to account for finite resolution, and because we computed the eigenstate spectrum with higher accuracy than the wavefunction propagation. Convergence tests show that the eigenstates for levels up to ∼1030 cm⁻¹ are converged to less than 1.5 cm⁻¹, much less than the errors of the PES representation (see ESI† for further details). This is the first time that so many states have been computed to such a high accuracy for such a high-dimensional, fluxional cluster.

Fig. 1 Infrared spectrum of the Zundel ion: the experimental predissociation spectrum using neon tagging (gray)^13,17 and computed spectra using either the BBSM (red) or the HBB (blue, negative intensities) potential energy surfaces; the inset magnified the spectral region in the gray box. The computed spectra are either based on time-dependent wavefunction propagations (td-TTNS: lines) or based on eigenstate optimizations (ti-TTNS: sticks). The wavefunction labels mark the states as analyzed in the text. All eigenstates computed, regardless of intensities, are shown as small lines on the abscissa.

While both PESs, in particular the newer BBSM one, yield spectra that overall are in very good agreement with the experimental spectrum, a closer look reveals a striking difference between the simulated spectra: the most characteristic signal of the IR spectrum, the doublet around 1000 cm⁻¹, consists of two dominant peaks on the BBSM PES whereas on the HBB PES it consists of three peaks! This is in contrast to the pioneering MCTDH simulations on the HBB PES.^14,15,18,32 These early MCTDH simulations of such a complicated system were less accurate and thus had simply missed the additional peak on the HBB PES. Indeed, our modern ML-MCTDH simulations with improved convergence parameters also yield three peaks on the HBB PES and overall are in full agreement with our TTNS simulations. Comparing the overlaps between the states on each PES reveals that the third peak observed for the HBB PES is not missing in the BBSM PES, but it is blueshifted and has a much smaller intensity.

We focus now on the nature of the three dominating vibrational transitions and pinpoint the differences resulting from the quality of the PESs. According to Fig. 1, we label these states as Ψ_a, Ψ_b, and Ψ_c, respectively, see also Fig. 1. When needed, we will use superscripts to display which PES has been used for their optimization. There are two additional, satellite peaks between Ψ^BBSM_b and Ψ^BBSM_c that have a similar intensity than Ψ^BBSM_c and contain a combination of various excitations. More details on these and other eigenstates that do not dominate the doublet will be presented elsewhere.

Our main findings, namely (i) the nature of the signals dominant in the doublet, (ii) the decomposition of the corresponding wavefunctions into coupling zero-order states, and (iii) their critical dependence on subtle energetic shifts, are summarised in Fig. 2. Panel (a) displays the coupling – a Fermi resonance – identified in the previous MCTDH studies on the HBB PES.^14,19 There, the doublet consists of two states that can be described as linear combination (entanglement) of two zero-order states: one zero-order state with one quantum in the water–water stretch motion (R) and two quanta in the two wagging modes (|02〉–|20〉), and another zero-order state with one quantum in the proton transfer motion (z). These two zero-order states are labeled |1R, 02–20〉 and |1z〉 respectively. Panel (b) displays our revised coupling scheme on the HBB PES. There, an additional zero-order state |04–40〉, consisting of four quanta in the wagging motion, has a low enough energy that it couples with the two other zero-order states, leading to a more complex quantum resonance and a significant intensity sharing of the |Ψ^HBB_b〉 and |Ψ^HBB_c〉. Panel (c) displays the coupling scheme on the BBSM PES, which, ironically, is much more similar to that of panel (a) than to (b). Here, the |04–40〉 state is blueshifted and, therefore, does not contribute anymore to the doublet. There is no intensity-sharing and |Ψ^BBSM_c〉 only appears in the spectrum as one of two satellite peaks on the blue wing of |Ψ^BBSM_b〉 – in good agreement with the experimental IR spectrum as depicted in the upper panel of Fig. 1. (The other satellite peak is more complicated and will be discussed elsewhere.)

Fig. 2 Schematics of the coupling scheme of the zero-order states around the doublet in the Zundel infrared spectrum. Panel (a) displays the previously identified coupling on the pioneering HBB PES.^14,19 Panel (b) displays our revised coupling on the same PES. Panel (c) displays the coupling on the recent BBSM PES. The insets visualize schematically the vibrational excitations of the zero-order states. The energies of the zero-order states are estimated either from the exact |1R〉 and |02–20〉/|04–40〉 states, or from the |1z, 1α〉 combination state, where α represents the torsion motion.

Let us now explain our findings in more detail. Our analysis is based on the time-independent eigenstate simulations, which enable us to directly analyze the wavefunctions responsible for each peak, thus providing rigorous assignment of IR peaks to the structural dynamics at the level of vibrational motion. Moreover, we confirm this analysis by computing overlaps of the wavefunctions with constructed zero-order states (see ESI†). Fig. 3 shows cuts of the three most relevant eigenstates Ψ_a, Ψ_b, and Ψ_c on the two PESs along one of the wagging motions and the proton transfer motion. Each node (zero-crossing) corresponds to one quantum in a particular coordinate. The plots of Ψ_a clearly indicate two quanta (nodes) for the wagging motion (|02–20〉) and one quantum for the proton transfer motion (|1z〉). Further inspections reveal an additional quantum for the water–water stretch motion (|1R〉). This confirms the assignment shown in Fig. 2: |Ψ_a〉 consists of two entangled zero-order states, namely |02–20; 1R〉 and |1z〉. There is no significant difference between Ψ^BBSM_a and Ψ^HBB_a.

Fig. 3 Representative cuts of the three relevant wavefunctions Ψ_a, Ψ_b and Ψ_c corresponding either to the doublet and a satellite peak (BBSM, upper panels), or to the triplet (HBB PES, lower panels) in the IR spectrum. The abscissa shows the wagging (pyramidalization) motion of one of the water molecules whereas the ordinate shows the proton transfer motion. The red lines denote the zero contours.

In stark contrast to Ψ_a, for Ψ_b and Ψ_c we find significant differences between the BBSM PES and the HBB PES. At first sight, Ψ^BBSM_b looks much more similar to Ψ^HBB_c than to Ψ^HBB_b. We first analyze the states on the BBSM PES and then compare them with the states on the HBB PES. Ψ^BBSM_b has a similar nodal pattern as Ψ_a, as it also displays two quanta along the wagging motion and one quantum along the proton transfer motion. These similarities are the characteristics of a Fermi resonance. Ψ^BBSM_c, however, shows a dominant contribution of not two, but four quanta in the wagging motion. The rich vibrational structure along the two wagging coordinates is displayed in Fig. 4. This enables us to clearly assign Ψ^BBSM_c as |04–40〉 wagging state.

Fig. 4 Representative cuts of the two wavefunctions Ψ_b and Ψ_c along the two wagging motions computed with the BBSM (upper panels) and with the HBB (lower panels) PESs. The red lines denote the zero contours.

We now turn to the states on the HBB PES. For Ψ^HBB_b and unlike Ψ^BBSM_b, in addition to one quantum in the proton transfer motion, there are not two but four quanta along the wagging motion. Fig. 4 clearly shows the similarities to Ψ^BBSM_c. Here, however, Fig. 3 and additional analysis reveals a more dominant excitation along the proton transfer motion, and we assign Ψ^HBB_b as linear combination of the |04–40〉 wagging and the |1z〉 proton transfer zero-order states. Finally, similar to Ψ_a, we identify in Ψ^HBB_c not only dominating contributions of the |02–20; 1R〉 and the |1z〉 states, but also contributions of the |04–40〉 wagging state. Thus, on the HBB PES next to the resonance of the |02–20; 1R〉 and the |1z〉 state, which alone would lead to two dominant peaks, the |1z〉 state also resonates with the |04–40〉 wagging state. This coupling creates an additional splitting, leading to overall three peaks in the spectrum, a clear example of intensity sharing.

Why do the three particular zero-order states, |1R, 02–20〉, |04–40〉, and |1z〉, couple on the HBB PES and why is the |04–40〉 not involved in the coupling on the BBSM PES? Estimating the energetics of the zero-order states as shown in Fig. 2 is the clue. The |02–20〉 and the |1R〉 states alone do not have enough energy (374 and 546 cm⁻¹, respectively) to couple with |1z〉. Only a combination of these vibrational excitations leads to a state with similar energy than the |1z〉 state, which is required for efficient coupling. Next to the |02–20; 1R〉 state, on the HBB PES the |04–40〉 state has an energy that is very close to that of |1z〉. This leads to additional coupling, resulting in Ψ^HBB_b and Ψ^HBB_c. While |04–40〉 and |02–20; 1R〉 alone are barely IR active, large components of the |1z〉 lead to intensity sharing and thus large IR intensities in the spectrum for all three entangled states.

In contrast, on the BBSM PES the |04–40〉 state is blueshifted by 40 cm⁻¹ and the |1z〉 state is redshifted by 115 cm⁻¹. This disfavors additional coupling between these states. Hence, only the |02–20; 1R〉 couples with |1z〉, leading to the pair of states, Ψ^BBSM_a and Ψ^BBSM_b, which form a classical Fermi resonance. Ψ^BBSM_c is almost a pure |04–40〉 state with only minor contributions from |1z〉 and thus a weak IR intensity.

Which of these two scenarios reflects now the actual situation in experiment? While the experimental resolution is not high enough to fully reject the possibility of two states lying under one of the peaks of the doublet, there are two indications that the scenario on the BBSM PES is more realistic: (1) the peak positions obtained from the newer BBSM PES are closer to the experimental spectrum and the PES is much closer to basis-set-extrapolated CCSD(T) energies as shown in the left panels of Fig. S3.† (2) Ψ^BBSM_c can be attributed to one of two satellite peaks in the experimental spectrum, and (3) another state between Ψ^BBSM_b and Ψ^BBSM_c likely corresponds to the second satellite peak. We thus conclude that the additional coupling seen in the older HBB PES is most likely an artifact that has previously been missed due to less accurate computation of the IR spectrum, which is only now disclosed in view of much improved methodologies.

Based on these findings, it is striking to demonstrate that such relatively small changes in the energetics lead to such drastic differences in the entanglement of the eigenstates and in the IR intensities, which serve as the experimental observables for these intricate phenomena. Moreover, we infer that minute changes in the environment can lead to even larger energy differences and we anticipate effects similar to those shown here for microsolvated clusters and larger clusters with additional solvation shells. To give three examples: (1) deuteration of the Zundel ion significantly alters the energetics of the zero-order states and vastly complicates the IR spectrum.¹⁸ (2) The experimental spectrum of H₅O₂⁺·Ar (and that of H₅O₂⁺·H₂) is strikingly different to that of H₅O₂⁺·Ne, as it displays four and not two peaks around 1000 cm⁻¹,¹⁷ indicating a similar complex coupling situation (the Ar atom attaches to one of the OH units and thus lifts degeneracy; in contrast, the spectra of H₅O₂⁺·Ne and H₅O₂⁺·He do not show significant differences and their spectra should be very close to that of the bare H₅O₂^+ 34). (3) Recent ML-MCTDH simulations on the solvated hydronium or Eigen ion (H₉O₄⁺) reveal the contribution of dozens of eigenstates that dominate the IR activation of the hydronium O–H stretch motion and show that the proton vibrations of the Eigen ion can be understood in terms of an embedded Zundel subunit.³⁵ Likewise, recent experimental studies indicate that the dynamics of protonated water clusters can be related to fluctuations of local electrical fields,³⁶ which also appear for solvated H₅O₂⁺.^9,37,38

4 Conclusions

In conclusion, we report high-accuracy, state-resolved simulations of the Zundel ion on two PESs with unprecedented accuracy. Using variational DMRG-like algorithms adapted to tree tensor network states and polyspherical coordinates, we are able to compute up to 900 eigenstates up to 1900 cm⁻¹ with hitherto unreached accuracy for such a complex. Obtaining this large number of eigenstates has so far never been accomplished for a 15-dimensional, fluxional system using curvilinear coordinates, and this will serve as a benchmark for future quantum dynamics method developments. Only applying this advanced methodology allows us to draw definite conclusions on the quantum dynamics of the Zundel ion, in particular for the most dominant feature of the IR spectrum, namely the doublet around 1000 cm⁻¹. Previous assignments are shown to be incomplete due to several limitations of the methods available close to two decades ago. The novel assignment of the doublet around 1000 cm⁻¹ not only confirms the Fermi resonance nature of this proton transfer mode, but also explains so-far unassigned low-intensity satellite peaks in the experimental spectrum. In broader terms, our findings not only highlight the striking mechanistic consequences of the strong anharmonicity and fluxional character of this challenging but important molecular cluster, but also reveal how subtle effects of the energetics of the zero-order states can lead to qualitative changes of the IR spectrum, revealing fundamental quantum effects such as resonances, entanglement and intensity sharing. We are therefore convinced that the same computational approach as unleashed here for the first time will unveil similar spectral complexity upon assigning high-resolution IR spectra in a multitude of molecular complexes, notably when it comes to proton transfer in (micro-)solvation environments.

Data availability

Further data is available from the authors upon reasonable request.

Author contributions

H. R. L. and O. V. conceived the study; H. R. L. designed the study, developed the TTNS methodology and performed and analyzed the TTNS simulations; M. S. and O. V. performed the MCCPD fits and additional ML-MCTDH simulations; R. B., F. B., C. S., and D. M. created the interface for the BBSM PES, provided samples for the MCCPD procedure, and performed benchmarks for DMC configurations provided by M. S.; H. R. L. drafted the manuscript; all authors contributed to discussing the results and editing the manuscript.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank Professor Joel Bowman for having made their potential energy surface available. HRL acknowledges support from the University of California Merced start-up funding during the last part of this work. HRL acknowledges support from a postdoctoral fellowship from the German Research Foundation (DFG) via grant LA 4442/1-1 during the first part of this work. HRL acknowledges computational time both on the Pinnacles cluster at UC Merced (supported by NSF OAC-2019144) and at the Resnick High Performance Computing Center, a facility supported by the Resnick Sustainability Institute at the California Institute of Technology. The Bochum work has been funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy – EXC 2033 – 390677874 – RESOLV as well as by the individual DFG grant MA 1547/19 to DM and supported by the “Center for Solvation Science ZEMOS” funded by the German Federal Ministry of Education and Research and by the Ministry of Culture and Research of North Rhine-Westphalia. RB acknowledges funding from the Studienstiftung des deutschen Volkes and CS acknowledges partial financial support from the Alexander von Humboldt-Stiftung. MS and OV thank the High Performance Computing Center in Stuttgart (HLRS) under the grant number HDQM_MCT as well as the bwHPC project of the state of Baden-Württemberg under grant number bw18K011 for providing computational resources.

Notes and references

1. K. R. Asmis, N. L. Pivonka, G. Santambrogio, M. Brümmer, C. Kaposta, D. M. Neumark and L. Wöste, Science, 2003, 299, 1375–1377.
2. J.-W. Shin, N. I. Hammer, E. G. Diken, M. A. Johnson, R. S. Walters, T. D. Jaeger, M. A. Duncan, R. A. Christie and K. D. Jordan, Science, 2004, 304, 1137–1140.
3. J. M. Headrick, E. G. Diken, R. S. Walters, N. I. Hammer, R. A. Christie, J. Cui, E. M. Myshakin, M. A. Duncan, M. A. Johnson and K. D. Jordan, Science, 2005, 308, 1765–1769.
4. J. A. Fournier, C. J. Johnson, C. T. Wolke, G. H. Weddle, A. B. Wolk and M. A. Johnson, Science, 2014, 344, 1009–1012.
5. C. T. Wolke, J. A. Fournier, L. C. Dzugan, M. R. Fagiani, T. T. Odbadrakh, H. Knorke, K. D. Jordan, A. B. McCoy, K. R. Asmis and M. A. Johnson, Science, 2016, 354, 1131–1135.
6. D. Marx, M. E. Tuckerman, J. Hutter and M. Parrinello, Nature, 1999, 397, 601–604.
7. H. J. Zeng and M. A. Johnson, Annu. Rev. Phys. Chem., 2021, 72, 667–691.
8. M. Thämer, L. De Marco, K. Ramasesha, A. Mandal and A. Tokmakoff, Science, 2015, 350, 78–82.
9. F. Dahms, B. P. Fingerhut, E. T. J. Nibbering, E. Pines and T. Elsaesser, Science, 2017, 357, 491–495.
10. J. A. Fournier, W. B. Carpenter, N. H. C. Lewis and A. Tokmakoff, Nat. Chem., 2018, 10, 932–937.
11. B. Dereka, Q. Yu, N. H. C. Lewis, W. B. Carpenter, J. M. Bowman and A. Tokmakoff, Science, 2021, 371, 160–164.
12. M. E. Tuckerman, D. Marx, M. L. Klein and M. Parrinello, Science, 1997, 275, 817–820.
13. N. I. Hammer, E. G. Diken, J. R. Roscioli, M. A. Johnson, E. M. Myshakin, K. D. Jordan, A. B. McCoy, X. Huang, J. M. Bowman and S. Carter, J. Chem. Phys., 2005, 122, 244301.
14. O. Vendrell, F. Gatti and H.-D. Meyer, Angew. Chem., Int. Ed., 2007, 46, 6918–6921.
15. O. Vendrell, F. Gatti and H.-D. Meyer, J. Chem. Phys., 2007, 127, 184303.
16. M. Baer, D. Marx and G. Mathias, Angew. Chem., Int. Ed., 2010, 49, 7346–7349.
17. L. R. McCunn, J. R. Roscioli, M. A. Johnson and A. B. McCoy, J. Phys. Chem. B, 2008, 112, 321–327.
18. O. Vendrell, F. Gatti and H.-D. Meyer, Angew. Chem., Int. Ed., 2009, 48, 352–355.
19. O. Vendrell, F. Gatti and H.-D. Meyer, J. Chem. Phys., 2009, 131, 034308.
20. S. M. Dietrick and S. S. Iyengar, J. Chem. Theory Comput., 2012, 8, 4876–4890.
21. G. Bertaina, G. Di Liberto and M. Ceotto, J. Chem. Phys., 2019, 151, 114307.
22. X. Huang, B. J. Braams and J. M. Bowman, J. Chem. Phys., 2005, 122, 044308.
23. C. Schran, J. Behler and D. Marx, J. Chem. Theory Comput., 2020, 16, 88–99.
24. R. Beckmann, F. Brieuc, C. Schran and D. Marx, J. Chem. Theory Comput., 2022 DOI:10.1021/acs.jctc.2c00511.
25. H. R. Larsson, J. Chem. Phys., 2019, 151, 204102.
26. M. Schröder and H.-D. Meyer, J. Chem. Phys., 2017, 147, 064105.
27. M. Schröder, J. Chem. Phys., 2020, 152, 024108.
28. S. R. White, Phys. Rev. Lett., 1992, 69, 2863–2866.
29. S. R. White, Phys. Rev. B, 1993, 48, 10345–10356.
30. J. Haegeman, C. Lubich, I. Oseledets, B. Vandereycken and F. Verstraete, Phys. Rev. B, 2016, 94, 165116.
31. F. A. Y. N. Schröder, D. H. P. Turban, A. J. Musser, N. D. M. Hine and A. W. Chin, Nat. Commun., 2019, 10, 1062.
32. O. Vendrell, M. Brill, F. Gatti, D. Lauvergnat and H.-D. Meyer, J. Chem. Phys., 2009, 130, 234305.
33. H. Wang and M. Thoss, J. Chem. Phys., 2003, 119, 1289–1299.
34. C. J. Johnson, A. B. Wolk, J. A. Fournier, E. N. Sullivan, G. H. Weddle and M. A. Johnson, J. Chem. Phys., 2014, 140, 221101.
35. M. Schröder, F. Gatti, D. Lauvergnat, H.-D. Meyer and O. Vendrell, arXiv:2204.03744, 2022.
36. N. Yang, C. H. Duong, P. J. Kelleher and M. A. Johnson, Nat. Chem., 2020, 12, 159–164.
37. F. Dahms, R. Costard, E. Pines, B. P. Fingerhut, E. T. J. Nibbering and T. Elsaesser, Angew. Chem., Int. Ed., 2016, 55, 10600–10605.
38. W. B. Carpenter, Q. Yu, J. H. Hack, B. Dereka, J. M. Bowman and A. Tokmakoff, J. Chem. Phys., 2020, 153, 124506.

---

Footnotes

† Electronic supplementary information (ESI) available. See https://doi.org/10.1039/d2sc03189b

‡ Current address: Laboratoire Matière en Conditions Extrêmes, Université Paris-Saclay, CEA, DAM, DIF, 91297 Arpajon, France.

§ Current address: Yusuf Hamied Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW, UK.

---