# Supplementary information:
Proton-transfer spectroscopy beyond the normal-mode scenario

Florian N. Brünig  Freie Universität Berlin, Department of Physics, 14195
Berlin, Germany  Paul Hillmann  Freie Universität Berlin, Department of
Physics, 14195 Berlin, Germany  Won Kyu Kim  Korea Institute for Advanced
Study, School of Computational Sciences, Seoul 02455, Republic of Korea  Jan
O. Daldrop  Freie Universität Berlin, Department of Physics, 14195 Berlin,
Germany  Roland R. Netz<EMAIL_ADDRESS>Freie Universität Berlin,
Department of Physics, 14195 Berlin, Germany

(August 28, 2024)

##  I Quantum zero-point motion effects

Quantum zero-point motion is known to smear out particle distributions and
thereby to increase the particle density at barriers. Therefore in this
section, we estimate the potential barrier height of a double-well potential
for which the ground-state probability density develops a single minimum.

The one-dimensional stationary Schrödinger equation

$\displaystyle-\frac{\hbar^{2}}{2m}\frac{\partial^{2}}{{\partial
d}^{2}}\phi(d)$ $\displaystyle=(U(d)-E)\phi(d)$ (S1)

is solved numerically for the quartic double-well potential

$\displaystyle U(d)$ $\displaystyle=U_{0}\left(\left(\frac{d}{d_{\rm
TW}^{*}}\right)^{2}-1\right)^{2},$ (S2)

which is shown in fig. S1A to approximate the effective potential describing
the excess proton distribution obtained from our ab initio molecular-dynamics
(AIMD) simulations very well.

\begin{overpic}[width=212.47617pt]{{/../figs3/fe/zundel_ffit_dw}.pdf}
\put(-2.0,56.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs3/fe/zundel_doverU}.pdf}
\put(-2.0,56.0){\huge\bf B} \end{overpic}

Figure S1: A Free-energy profiles $U(d)$ for fixed $R_{\mathrm{OO}}$,
extracted from constrained ab initio molecular-dynamics (AIMD) simulations of
the H5O${}_{2}{}^{+}$ cation, are shown as colored lines. Fits according to
eq. (S2) are shown as black dotted lines with fit parameters given in the
legend. B Dependence of $d^{*}_{\rm TW}$ on $U_{0}$. Here a linear fit is
shown as a black dashed line.

The equations are rescaled as

$\displaystyle\frac{\partial^{2}}{{\partial\tilde{x}}^{2}}\phi(\tilde{x})$
$\displaystyle=(\widetilde{U}(\tilde{x})-\widetilde{E})\phi(\tilde{x}),$ (S3)
$\displaystyle\widetilde{V}(\tilde{x})$
$\displaystyle=\widetilde{U}_{0}(\tilde{x}^{2}-1)^{2},$ (S4)

with $\tilde{x}=d/d_{\rm TW}^{*}$, $\widetilde{U}_{0}=2mU_{0}{d_{\rm
TW}^{*}}^{2}/\hbar^{2}$ and $\widetilde{E}=2mE{d_{\rm TW}^{*}}^{2}/\hbar^{2}$.
Using the proton mass
$m=$1.7\text{\times}{10}^{-27}\text{\,}\mathrm{k}\mathrm{g}$$,
$\hbar=$1.1\text{\times}{10}^{-34}\text{\,}\mathrm{J}\mathrm{s}$$ and
$d^{*}_{\rm TW}=$0.22\text{\,}\mathrm{\mbox{\AA}}$$ for the system with
$R_{\mathrm{OO}}=$2.64\text{\,}\mathrm{\mbox{\AA}}$$, the rescaling factor for
the energies is $\hbar^{2}/(2m{d_{\rm
TW}^{*}}^{2})=$7.4\text{\times}{10}^{-21}\text{\,}\mathrm{J}$=1.8\ k_{B}T$ at
$300\text{\,}\mathrm{K}$. Note that this scaling factor implicitly depends on
$U_{0}$ due to the dependence of $d^{*}_{\rm TW}$ on $U_{0}$ shown in fig.
S1B.

\begin{overpic}[width=212.47617pt]{{/../figs/groundstate_energy_over_barrier_alt_symbol2}.png}
\put(-2.0,56.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs/groundstate_curv_over_barrier}.png}
\put(-2.0,56.0){\huge\bf B} \end{overpic}

Figure S2: A Dimensionless ground state energy $\widetilde{E}$ as a function
of the barrier height $\widetilde{U}_{0}$ (green dots). A linear fit according
to $\widetilde{E}(U_{0})=0.601\,\widetilde{U}_{0}+0.542$ is shown as a black
dashed line. The crossing between the ground-state energy $\widetilde{E}$ and
$\widetilde{E}=\widetilde{U}_{0}$ (red line) is marked as a vertical blue line
at
$\widetilde{U}_{c}:=\widetilde{U}_{0}(\widetilde{E}=\widetilde{U}_{0})=1.36$.
B Curvature of the ground state probability density at the origin as a
function of $\widetilde{U}_{0}$. A linear fit according to
$\partial^{2}\phi^{2}/{\partial\tilde{x}}^{2}(\tilde{x}=0)=0.334\,\widetilde{U}_{0}-0.455$
is shown as a black dashed line. The point where the curvature vanishes is
marked as a vertical blue line, which is given by $\widetilde{U}_{c}=1.36$.

Dimensionless ground state energies $\widetilde{E}$ are numerically computed
for $\widetilde{U}_{0}\in[1,2]$ using bisection to find the eigenvalues of the
dimensionless Schrödinger equation eq. (S3). The ground state energy
$\widetilde{E}$ shows approximately linear dependence on
$\widetilde{U}_{0}\in[1,2]$ (see green dots in fig. S2A). A minimum of the
ground-state density is defined by a positive curvature at the origin,
$\frac{\partial^{2}}{{\partial x}^{2}}\phi^{2}(\tilde{x}=0)>0$. From the
Schrödinger equation eq. (S3) it becomes clear that the wave function of the
ground state develops a minimum at the origin for
$\widetilde{E}<\widetilde{U}_{0}$. It follows for the curvature of the
probability density

$\displaystyle\frac{\partial^{2}}{{\partial\tilde{x}}^{2}}\phi^{2}$
$\displaystyle=2\left[\left(\frac{\partial\phi}{\partial\tilde{x}}\right)^{2}+(\widetilde{U}(\tilde{x})-\widetilde{E})\phi^{2}\right].$
(S5)

Since the ground state is symmetric around the origin this becomes at the
origin

$\frac{\partial^{2}}{{\partial\tilde{x}}^{2}}\phi^{2}=2(\widetilde{U}_{0}-\widetilde{E})\phi^{2}.$
(S6)

Fig. S2B shows the curvature at the origin $\frac{\partial^{2}}{{\partial
x}^{2}}\phi^{2}(\tilde{x}=0)$ as a function of $\widetilde{U}_{0}$ (green
dots). A linear fit is used to determine the critical potential strength at
which the curvature vanishes as $\widetilde{U}_{c}=1.36$. Fig. S3 shows a
comparison of the ground state density distributions for
$\widetilde{U}_{0}=\\{1,\widetilde{U}_{c},1.8\\}$.

From this analysis we conclude that when treating the excess proton quantum-
mechanically, a minimum in the ground state probability density, reflecting
the effect of the barrier, appears for barrier heights of the double-well
potential

$\displaystyle U_{0}>\frac{\hbar^{2}}{2m{{d_{\rm
TW}^{*}}^{2}(U_{0})}}\,\widetilde{U}_{c}\,k_{B}T,$ (S7)

which holds for $U_{0}>2.5\,k_{B}T$ as can be read off fig. S4. Therefore,
using the naive assumption that a minimum in the probability density at the
barrier top indicates the presence of an effective quantum-mechanical barrier,
which is far from obvious, one could speculate that the presented results for
the spectroscopic signatures of proton barrier-crossing are expected to
survive quantum-mechanical zero-point motion effects for high enough potential
barriers.

\begin{overpic}[width=212.47617pt]{{/../figs/groundstate_wf}.png}
\put(-2.0,56.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs/groundstate_wf_zoom}.png}
\put(-2.0,56.0){\huge\bf B} \end{overpic}

Figure S3: A Ground state probability densities for $\widetilde{U}_{0}=1$,
$\widetilde{U}_{0}=\widetilde{U}_{c}$ and $\widetilde{U}_{0}=1.8$. B Ground
state probability densities around the origin.
\begin{overpic}[width=212.47617pt]{{/../figs3/fe/zundel_UcoverU}.pdf}
\put(-2.0,56.0){} \end{overpic} Figure S4: A minimum in the ground state
probability density appears for $\frac{2m}{\hbar^{2}k_{B}T}{{d_{\rm
TW}^{*}}^{2}(U_{0})}U_{0}>\widetilde{U}_{c}$ as a rearrangement of eq. (S7)
with $U_{0}$ being the barrier height of the double-well potential,
$\widetilde{U}_{c}=1.36$ taken from fig. S2B and ${d_{\rm TW}^{*}(U_{0})}$
taken from the linear fit in fig. S1B.

##  II Summary of previous experimental and theoretical studies on infrared
spectra of the H5O${}_{2}{}^{+}$ cation

\begin{overpic}[width=212.47617pt]{{/../figs/zundel_compare_exp}.pdf}
\put(-4.0,62.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs/zundel_compare_theo}.pdf}
\put(-4.0,95.0){\huge\bf B} \end{overpic}

Figure S5: Collection of experimental (A) and theoretical (B) IR spectra of
the H5O${}_{2}{}^{+}$ cation in the proton-transfer regime (colored lines),
compared to the IR spectrum obtained from AIMD simulations in this study,
shown in black. A The experimental spectra were recorded using multiple-
photon-dissociation (IRPMD) spectroscopy [1, 2], predissociation (IRPD)
spectroscopy in Argon [3] and Neon [4] and Fourier-transform IR (FTIR)
spectroscopy of H5O${}_{2}{}^{+}$ cations solvated in acetonitrile [5]. B
Theoretical spectra were obtained using AIMD simulations on the MP2/cc-pVTZ
level [6], the multiconfiguration time-dependent Hartree (MCTDH) method [7],
CPMD simulations using the BLYP functional [8], CPMD simulations using the
BLYP functional [9], MD simulations on the CCSD(T) potential energy surface
using MP2 dipole moment functions (based on [10]) [9], AIMD simulations on the
BLYP-D3 TZV2P level [11], MD simulations using the multistate-empirical-
valence-bond (MS-EVB) method [12], AIMD simulations on the DZVP-BLYP level
[13], classical MD, centroid MD (CMD), as well as ring polymer MD simulations
(RPMD, TRPMD)[14] on the CCSD(T) potential energy surface [10] and normal mode
analysis on the 6-311++G(d,p)-B3LYP level performed on states obtained from a
multistate-empirical-valence-bond (MS-EVB) MD simulation [15]. All theoretical
results were obtained for a single H5O${}_{2}{}^{+}$ cation, except for [15],
who computed the spectra from clusters of 16–18 water molecules and an excess-
proton.

A collection of IR spectra of the H5O${}_{2}{}^{+}$ cation in the so-called
proton-transfer regime,
$600\text{\,}\mathrm{c}\mathrm{m}^{-1}1500\text{\,}\mathrm{c}\mathrm{m}^{-1}$,
obtained from experiments is given in fig. S5A and from theoretical
calculations in fig. S5B. The spectra are compared to the spectrum of the
unconstrained H5O${}_{2}{}^{+}$ cation obtained from AIMD simulations in this
study (black lines on top). The spectra vary greatly among different
experimental studies in this regime, highlighting the subtle differences
between experimental techniques and also pointing to a pronounced temperature
dependence of the spectra. Fridgen et al. [2] say that “at this point, we are
unable to determine the source of the discrepancy between the present infrared
spectrum and that obtained by Asmis et al. [1]”. Sauer and Döbler [6] discuss
“that differences between the IRMPD spectrum [1] and the IRPD spectrum [3]
should be due to the different excitation mechanisms and/or different
temperatures”. Hammer et al. [4] state “the IR profiles obtained in these two
measurements [1, 2] were markedly different, perhaps reflecting the different
ion sources used in the two experiments and/or the specific fluence
characteristics of the laser sources”. Vendrell et al. [7] introduce the topic
stating that “spectra could not be consistently assigned in terms of
fundamental frequencies and overtones of harmonic vibrational modes due to
large amplitude anharmonic displacements and couplings of the cluster [1, 2,
16, 4]” and conclude their study by pointing out that their “reported
calculations are in excellent agreement to the experimental measurements of
Refs. [16, 4] on the predissociation spectrum of H5O2+·Ne.”, which is the most
widely accepted low-temperature spectrum to date. Furthermore, nuclear quantum
effects were studied in detail by Rossi et al. [14], who compared spectra from
different ring polymer MD simulations, and Park et al. [8] and Baer et al.
[11] successfully investigated the messenger-induced changes of spectra
apparent in IRMPD techniques. The temperature dependence was analyzed by Park
et al. [8] and Kaledin et al. [9] in Car-Parinello MD (CPMD) simulations,
presented in fig. S6, and classical MD simulations on the CCSD(T) potential
energy surface using MP2 dipole moment functions (based on [10]), stating that
“classical MD simulations at a temperature of $30\text{\,}\mathrm{K}$
qualitatively reproduce many of the key features in the experimental
vibrational [Ar-] predissociation”. More recent studies have focused on the
H5O${}_{2}{}^{+}$ cation in larger water clusters or in bulk water or other
solvents and at room temperature [17, 13, 5, 15, 12, 18]. These studies
identify remarkably broad IR features in the
$1000\text{\,}\mathrm{c}\mathrm{m}^{-1}1200\text{\,}\mathrm{c}\mathrm{m}^{-1}$
regime, associated with “proton shuttling” [17], the “proton-transfer mode
(PTM)” [13, 15, 5] or the “shared proton stretch” [18].

The review of the published experimental and theoretical data shows that the
data produced in this study is in qualitative agreement with previous results.
However, especially in the region of
$800\text{\,}\mathrm{c}\mathrm{m}^{-1}1200\text{\,}\mathrm{c}\mathrm{m}^{-1}$
associated with the proton-transfer motion, large temperature and system-
dependent differences between the previously published results are observed
and reflect the incomplete understanding of proton-transfer dynamics.

\begin{overpic}[width=212.47617pt]{{/../figs/zundel_compare_theo_T}.pdf}
\end{overpic} Figure S6: Comparison of theoretical IR spectra of the
H5O${}_{2}{}^{+}$ cation in the proton-transfer regime at various temperatures
obtained from Car-Parinello MD (CPMD) simulations using the BLYP density
functional method [8, 9].

##  III Normal-mode analysis of the H5O${}_{2}{}^{+}$ cation

\begin{overpic}[width=433.62pt]{{/../figs/zundel_compare_spectra_nm_xyz_all}.pdf}
\put(2.0,79.0){\large\bf A} \put(37.0,79.0){\large\bf B}
\put(68.0,79.0){\large\bf C} \put(2.0,60.0){\large\bf D}
\put(37.0,60.0){\large\bf E} \put(68.0,60.0){\large\bf F}
\put(2.0,42.0){\large\bf G} \put(37.0,42.0){\large\bf H}
\put(68.0,42.0){\large\bf I} \put(2.0,24.0){\large\bf J}
\put(37.0,24.0){\large\bf K} \put(68.0,24.0){\large\bf L} \end{overpic} Figure
S7: IR spectra of the excess proton (colored lines) in the H5O${}_{2}{}^{+}$
cation for different constraint oxygen distances $R_{\mathrm{OO}}$ compared to
normal-mode spectra calculated from energetically optimized structures and
projected on the proton coordinate(grey dashed sticks). The discrete normal
modes sum to smooth spectra (grey lines) by including line-broadening at
finite temperature obtained from damped harmonic oscillations of the normal
modes according to eq. S36 in section X and assuming a friction coefficient of
$\gamma=$$16\text{\,}\mathrm{u}\mathrm{/}\mathrm{p}\mathrm{s}$. The spectra
are shown for different orientations (A–C, G–I: $x$-axis connecting the two
oxygens, D–F, J–L: $yz$-plane).

The normal modes are obtained as the Eigenvectors of the Hessian of the
optimal structure in a chosen coordinate system, the Eigenvalues give the
frequencies. The normal-mode analysis of the H5O${}_{2}{}^{+}$ systems was
performed for energetically optimized structures using the same set of
parameters as for the AIMD simulation.

A projection of the Eigenvectors on the proton coordinate allows to obtain the
relative magnitudes of the excess-proton spectra predicted from the normal-
mode analysis. These are scaled to the magnitudes of the excess-proton spectra
obtained directly from the AIMD simulations for various $R_{\mathrm{OO}}$ and
orientations in fig. S7. To improve the normal-mode spectra, line-broadening
of the normal modes at finite temperature is modeled by damped harmonic
oscillations (see SI section X). In the $yz$-plane the agreement is very good
for all barrier heights indicating, that the method is well suited for
modeling the excess-proton dynamics without a barrier. Along $x$ the normal
mode spectra serve to explain the major peaks of the spectra at high
frequencies which will be identified as the normal-mode contributions. The
broad low-frequency shoulder, which is related to the barrier-crossing
effects, can obviously not be modeled using normal modes.

The system with $R_{\mathrm{OO}}=$2.40\text{\,}\mathrm{\mbox{\AA}}$$
represents the global unconstrained minimum of the free-energy landscape,
shown in fig. LABEL:zundel_introD in the main text, and can therefore be used
to compare qualitatively to published data on the unconstrained
H5O${}_{2}{}^{+}$ cation, as done in tab. S1 for the range
$100\text{\,}\mathrm{c}\mathrm{m}^{-1}2000\text{\,}\mathrm{c}\mathrm{m}^{-1}$.
The normal-modes can be grouped to rocking and wagging modes in the range
$300\text{\,}\mathrm{c}\mathrm{m}^{-1}600\text{\,}\mathrm{c}\mathrm{m}^{-1}$,
the water-water stretch at around
$550\text{\,}\mathrm{c}\mathrm{m}^{-1}650\text{\,}\mathrm{c}\mathrm{m}^{-1}$,
which is missing in the constrained case, the widely discussed proton-transfer
mode along $x$, which varies greatly between
$800\text{\,}\mathrm{c}\mathrm{m}^{-1}1200\text{\,}\mathrm{c}\mathrm{m}^{-1}$,
a pair of modes between
$1300\text{\,}\mathrm{c}\mathrm{m}^{-1}1600\text{\,}\mathrm{c}\mathrm{m}^{-1}$,
associated with perpendicular excess-proton motion in the $yz$-plane, and the
two water-bending modes at
$1600\text{\,}\mathrm{c}\mathrm{m}^{-1}1800\text{\,}\mathrm{c}\mathrm{m}^{-1}$,
for which the lower one shows excess-proton motion in the $yz$-plane (named
in-phase or gerade) and the higher one shows excess-proton motion along $x$
(named out-of-phase or ungerade).

Illustrations of the normal modes in the range
$100\text{\,}\mathrm{c}\mathrm{m}^{-1}2000\text{\,}\mathrm{c}\mathrm{m}^{-1}$
of H5O${}_{2}{}^{+}$ cations with various contrained $R_{\mathrm{OO}}$ are
shown in fig. S8, where the colored frames indicate similar normal modes for
different values of $R_{\mathrm{OO}}$. Fig. S9 summarizes the frequencies of
the normal modes for all values of $R_{\mathrm{OO}}$ and highlights the
associated relative excess-proton-motion intensities either along $x$ or in
the $yz$-plane. Most modes shift only weakly, including the wagging and
rocking modes in the range
$300\text{\,}\mathrm{c}\mathrm{m}^{-1}600\text{\,}\mathrm{c}\mathrm{m}^{-1}$
and the pair of modes between
$1300\text{\,}\mathrm{c}\mathrm{m}^{-1}1600\text{\,}\mathrm{c}\mathrm{m}^{-1}$,
associated with perpendicular excess-proton motion in the $yz$-plane. However,
the proton-transfer mode ($1157\text{\,}\mathrm{c}\mathrm{m}^{-1}$ for
$R_{\mathrm{OO}}=$2.40\text{\,}\mathrm{\mbox{\AA}}$$) strongly increases for
values of $R_{\mathrm{OO}}\geq$2.40\text{\,}\mathrm{\mbox{\AA}}$$, as shown
previously by Wolke et al. [19] for the D5O${}_{2}{}^{+}$ cation, which
possibly explains why the proton-transfer mode in bulk is suggested to reside
at much higher wavenumbers compared to the gas-phase spectra [13]. The in-
phase water bending mode shows no shifting but the out-of-phase water bending
mode, associated with excess-proton motion along $x$, strongly shifts.
Generally, modes associated with proton motion along $x$ are shown to be
highly sensitive on the value of $R_{\mathrm{OO}}$, while modes that are
associated with proton motion in the $yz$-plane are not.

Table S1: Collection of previous identifications of normal modes in the H5O${}_{2}{}^{+}$ cation in the range $100\text{\,}\mathrm{c}\mathrm{m}^{-1}2000\text{\,}\mathrm{c}\mathrm{m}^{-1}$, compared to data obtained in this study with constrained $R_{\mathrm{OO}}=$2.40\text{\,}\mathrm{\mbox{\AA}}$$ in the first row. Bold fonts indicate significant motion of the excess proton (if reported). BLYP-D3, TZV2P, $R_{\mathrm{OO}}=$2.40\text{\,}\mathrm{\mbox{\AA}}$$ | wagging, rocking |  320, 426, 456, 470 |  |  |  p. trans. ($x$) | 1157 |  proton perp. (${yz}$) | 1350, 1428 | water bend | 1591, 1653
---|---|---|---|---|---|---|---|---|---|---
MP2, aug-cc-pVTZ [6] |  H2O wag, H2O rock | 362, 461, 532, 535 | sym. OHO str. | 627 | asym. OHO str. | 884 |  OHO bend $(y,z)$ | 1484, 1557 | HOH bend | 1710, 1765
MP2, cc-pVTZ [6] |  | 367, 451, 522, 529 |  | 631 |  | 961 |  | 1493, 1572 | (in/out phase) | 1710, 1773
MP2, cc-pVTZ(aug-O) [6] |  | 347, 456, 525, 532, |  | 626 |  | 904 |  | 1473, 1551 |  | 1706, 1761
CCSD(T), cc-pVTZ [6] |  | 350, 457, 531, 535 |  | 633 |  | 896 |  | 1493, 1572 |  | 1724, 1780
CCSD(T), MP2, aug-cc-pVTZ [10] |  | 339, 471, 532, 554 |  | 630 |  | 861 |  | 1494, 1574 |  | 1720, 1770
MP2, aug-cc-pVTZ [4] |  |  |  |  | OHO stretch | 808 |  |  | HOH bend | 1662, 1717
MCTDH [7] | wagging | 106, 108, 232, 254, 374, 422 | wat-wat str. | 550 | p. trans. doub. | 918, 1033 | proton perp. | 1391 | water bend | 1606, 1741
| rocking | 481, 915, 930, 943 |  |  |  |  | p. trans. + wat. str. | 1411 | (gerade, unger.) |
BLYP, aug-cc-pVTZ [8] |  |  |  |  | asym. OHO str. | 938 | OHO bend (y, z) | 1408, 1478 | HOH bend | 1647, 1702
B3LYP, aug-cc-pVTZ [8] |  |  |  |  |  | 891 |  | 1397, 1455 |  | 1619, 1678
MP2, aug-cc-pVTZ [8] |  |  |  |  |  | 868 |  | 1405, 1478 |  | 1627, 1680
BLYP-D3, TZV2P, AIMD [11] | wagging, rocking | 370, 431, 498, 514 | OO stretch | 571 |  p. trans. ($v_{q}$) | 971 |  protonyz ($\gamma_{p},\gamma_{p}^{\prime}$) | 1396, 1477 |  water bend ($\delta_{s},\delta_{a}$) |  1636, 1700
MS-EVB, EMA [12] |  |  |  |  |  | 1032 |  | 1311, 1560 |  | 1638, 1769
B3LYP, aug-cc-pVTZ [13] |  | 358 |  |  |  | 926 |  | 1446 |  | 1742

\begin{overpic}[width=346.89731pt]{{/../figs3/nm/zundel_nm_vis_d2.32}.pdf}
\put(0.0,22.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=346.89731pt]{{/../figs3/nm/zundel_nm_vis_d2.40}.pdf}
\put(0.0,22.0){\huge\bf B}
\end{overpic}\begin{overpic}[width=346.89731pt]{{/../figs3/nm/zundel_nm_vis_d2.48}.pdf}
\put(0.0,22.0){\huge\bf C}
\end{overpic}\begin{overpic}[width=346.89731pt]{{/../figs3/nm/zundel_nm_vis_d2.56}.pdf}
\put(0.0,22.0){\huge\bf D}
\end{overpic}\begin{overpic}[width=346.89731pt]{{/../figs3/nm/zundel_nm_vis_d2.64}.pdf}
\put(0.0,22.0){\huge\bf E}
\end{overpic}\begin{overpic}[width=346.89731pt]{{/../figs3/nm/zundel_nm_vis_d2.72}.pdf}
\put(0.0,22.0){\huge\bf F} \end{overpic}

Figure S8: Illustrations of the normal modes of the H5O${}_{2}{}^{+}$ cation
for various constrained $R_{\mathrm{OO}}$, projected on the $xy$ and $xz$
planes. The colored frames guide the eye through the shifting of a respective
normal mode through A-F.

\begin{overpic}[width=212.47617pt]{{/../figs3/nm/zundel_nm_allD_2_x}.pdf}
\put(0.0,62.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs3/nm/zundel_nm_allD_2_yz}.pdf}
\put(0.0,62.0){\huge\bf B} \end{overpic}

Figure S9: Normal mode frequencies of the excess-proton dynamics in the
H5O${}_{2}{}^{+}$ cation for various constrained $R_{\mathrm{OO}}$ as shown in
fig. S7. The shading indicates the relative amplitudes of the respective
normal modes, normalized to the maximal amplitude for each value of
$R_{\mathrm{OO}}$. The colors and grey dashed lines guide the eye through the
shifting of a respective normal mode, as also shown in fig. S8.

##  IV Decomposition of the H5O${}_{2}{}^{+}$ IR spectra into excess-proton
and water contributions

\begin{overpic}[width=138.76157pt]{{/../figs3/decomp/zundel_decomp_d2.32_f6_loglog}.pdf}
\put(-2.0,85.0){\large\bf A} \put(-2.0,45.0){\large\bf D}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/decomp/zundel_decomp_d2.40_f6_loglog}.pdf}
\put(-2.0,85.0){\large\bf B} \put(-2.0,45.0){\large\bf E}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/decomp/zundel_decomp_d2.48_f6_loglog}.pdf}
\put(-2.0,85.0){\large\bf C} \put(-2.0,45.0){\large\bf F}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/decomp/zundel_decomp_d2.56_f6_loglog}.pdf}
\put(-2.0,85.0){\large\bf G} \put(-2.0,45.0){\large\bf J}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/decomp/zundel_decomp_d2.64_f6_loglog}.pdf}
\put(-2.0,85.0){\large\bf H} \put(-2.0,45.0){\large\bf K}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/decomp/zundel_decomp_d2.72_f6_loglog}.pdf}
\put(-2.0,85.0){\large\bf I} \put(-2.0,45.0){\large\bf L} \end{overpic}

Figure S10: Decomposition of the IR spectra of the H5O${}_{2}{}^{+}$ cation
with two constrained oxygens at a given distance $R_{\mathrm{OO}}$ in
different directions (A–C, G–I: along the $x$-axis connecting the two oxygens,
D–F, J–L: along the $yz$-plane) into water and excess-proton contributions.
The total IR spectra are shown in grey, the spectra of the excess proton in
green and the spectra of the two water molecules in blue. The cross-
correlation spectrum defined as
$\widetilde{\chi}^{\prime\prime}_{\mathrm{H}^{+},\mathrm{H}_{2}\mathrm{O}}=\widetilde{\chi}^{\prime\prime}_{\mathrm{tot}}-\widetilde{\chi}^{\prime\prime}_{\mathrm{H}^{+}}-\widetilde{\chi}^{\prime\prime}_{\mathrm{H}_{2}\mathrm{O}}$
is shown in red. The dashed red line denotes negative values of the cross-
correlation spectrum.

By using Wannier centers for charge localization, the total dipole moment of
the simulation systems can be exactly decomposed into proton and water
contributions
$\mathbf{p}_{\mathrm{tot}}(t)=\mathbf{p}_{\mathrm{H}^{+}}(t)+\mathbf{p}_{\mathrm{H}_{2}\mathrm{O}}(t)$.
A comparison of the IR spectra of the total dipole moment
$\omega\widetilde{\chi}^{\prime\prime}_{\mathrm{tot}}$ of the
H5O${}_{2}{}^{+}$ cation for various constrained oxygen positions to the power
spectra of only the excess proton
$\omega\widetilde{\chi}^{\prime\prime}_{\mathrm{H}^{+}}$ and the power spectra
of the dipole moments of the two flanking water molecules
$\omega\widetilde{\chi}^{\prime\prime}_{\mathrm{H}_{2}\mathrm{O}}$ is shown in
fig. S10A-L in different directions. As discussed in SI section VII, the
cross-correlation spectra
$\omega\widetilde{\chi}^{\prime\prime}_{\mathrm{H}^{+},\mathrm{H}_{2}\mathrm{O}}$,
shown in red in fig. S10A-L, are proportional to the cross-correlations of
water-dipole-moment and excess-proton dynamics. Along $x$, the cross-
correlation spectra of the six systems shown here, are nearly entirely
positive as well as nearly proportional to the power spectrum of the proton
itself, indicating constructive coupling of the proton motion to the water
dipole moments along this axis, as previously shown [6]. An apparent exception
to this is the rocking and wagging regime at
$300\text{\,}\mathrm{c}\mathrm{m}^{-1}600\text{\,}\mathrm{c}\mathrm{m}^{-1}$
and the split OH-stretching mode at around
$3500\text{\,}\mathrm{c}\mathrm{m}^{-1}$, which shows a strong negative cross-
correlation spectrum in both $x$ direction and $yz$-plane. Adjacent to this
band all systems show two weak OH-stretching vibrations along $x$, likely a
Stark effect of the mean excess-proton field on the water dipole moments.

In the $yz$-plane the cross-correlation spectra are nearly entirely negative,
indicating mainly out-of-phase motion of the excess-proton and the water
dipole moments, producing generally weaker spectra in $yz$ compared to $x$.

It can be concluded that the IR spectrum of the H5O${}_{2}{}^{+}$ cation along
the axis connecting the two water oxygens reflects the excess-proton dynamics.
This justifies the IR signal to be used as a reporter of the excess-proton
dynamics, as well as the focus on its dynamics.

##  V Comparison of constrained and unconstrained dichroic IR spectra of the
H5O${}_{2}{}^{+}$ cation

\begin{overpic}[width=433.62pt]{{/../figs/zundel_compare_spectra_Boltzmann_xyz_all_x_2}.pdf}
\put(2.0,79.0){\large\bf A} \put(36.0,79.0){\large\bf B}
\put(67.0,79.0){\large\bf C} \put(2.0,60.0){\large\bf D}
\put(36.0,60.0){\large\bf E} \put(67.0,60.0){\large\bf F}
\put(2.0,42.0){\large\bf G} \put(36.0,42.0){\large\bf H}
\put(67.0,42.0){\large\bf I} \put(2.0,24.0){\large\bf J}
\put(36.0,24.0){\large\bf K} \put(67.0,24.0){\large\bf L} \end{overpic} Figure
S11: IR spectra of the H5O${}_{2}{}^{+}$ cation in different directions (A, D,
G, J: along the $x$-axis connecting the two oxygens, B, E, H, K: along the
$yz$-plane C, F, I, L: isotropic spectrum). A–I The colors correspond to
systems with different fixed oxygen-oxygen separation $R_{\mathrm{OO}}$, grey
lines denote the unconstrained system. J–L Comparison of the IR spectra of the
unconstrained system (grey lines) to a Boltzmann average (red lines) according
to eq. (S8) of the spectra in A-I. The HOH-bending mode and OH-stretching mode
of the water molecules are indicated in F and I. The short straight grey lines
indicate the $\omega^{2}$-scaling associated with the low-frequency shoulder
of the power spectra.

IR spectra of the unconstrained H5O${}_{2}{}^{+}$ cation and various systems
with constrained $R_{\mathrm{OO}}$ are shown in fig. S11 in different
directions. The $x$-axis corresponds to the direction connecting the two
oxygens, while the $yz$-plane is orthogonal to that axis, as illustrated in
fig. LABEL:zundel_introA in the main text. Interestingly, the spectra in
$x$-direction strongly depend on the value of the constrained oxygen
separation $R_{\mathrm{OO}}$ for frequencies lower than
$2500\text{\,}\mathrm{c}\mathrm{m}^{-1}$ (figs. S11A, S11D, S11G and S11J).
The HOH-bending mode is shifted to almost
$2000\text{\,}\mathrm{c}\mathrm{m}^{-1}$ for
$R_{\mathrm{OO}}=$2.32\text{\,}\mathrm{\mbox{\AA}}$$ (fig. S11A) and to
slightly lower frequencies than the unconstrained system for
$R_{\mathrm{OO}}\geq$2.40\text{\,}\mathrm{\mbox{\AA}}$$ (figs. S11D and S11G),
whereas the frequency of the OH-stretching mode is not affected by fixing
$R_{\mathrm{OO}}$.

The various spectra for the $yz$-plane in figs. S11B, S11E, S11H and S11K are
indistinguishable on the other hand. The dominant features are the librations
of the water molecules at $300\text{\,}\mathrm{c}\mathrm{m}^{-1}$ and the OH-
stretching mode at around $3500\text{\,}\mathrm{c}\mathrm{m}^{-1}$. While the
OH-stretching mode appears in both directions, the HOH-bending mode is mainly
visible along the $x$-axis. Note that in the $yz$-plane the HOH-bending mode
contributions of the outer hydrogens of each water molecule cancel out due to
symmetric motion with respect to the $x$-axis.

The isotropic spectra in figs. S11C, S11F, S11I and S11L depend on
$R_{\mathrm{OO}}$ only in the regime
$400\text{\,}\mathrm{c}\mathrm{m}^{-1}1700\text{\,}\mathrm{c}\mathrm{m}^{-1}$
since for lower frequencies, between
$10\text{\,}\mathrm{c}\mathrm{m}^{-1}400\text{\,}\mathrm{c}\mathrm{m}^{-1}$,
the $yz$-contributions to the isotropic spectrum are dominant. In particular,
this motivates the analysis of dynamics along $x$-direction, which shows the
dominant contribution to the isotropic spectrum except for very low
frequencies.

Interestingly, the spectra of the unconstrained systems are well recovered by
a Boltzmann average of the spectra of the systems with fixed $R_{\mathrm{OO}}$
according to

$\displaystyle\omega\widetilde{\chi}^{\prime\prime}(\omega)_{<R_{\mathrm{OO}}>}=\frac{\sum_{i}\omega\widetilde{\chi}^{\prime\prime}(\omega)_{R_{\mathrm{OO}_{i}}}e^{-U(R_{\mathrm{OO}_{i}})/k_{B}T}}{\sum_{i}e^{-U(R_{\mathrm{OO}_{i}})/k_{B}T}},$
(S8)

using the free energy along $R_{\mathrm{OO}}$ of the unconstrained system,
$U(R_{\mathrm{OO}})$, shown in fig. LABEL:zundel_introD in the main text. As
shown in figs. S11J, fig. S11K and fig. S11L, the agreement is very good along
all directions, which indicates a sufficient dynamic decoupling of the slow
oxygen coordinate $R_{\mathrm{OO}}$ from the proton coordinate $d$ with
respect to IR spectra and allows observations for the constrained systems to
be generalized to the unconstrained system.

##  VI Recrossing transfer paths of the excess-proton in the H5O${}_{2}{}^{+}$
cation

After crossing the barrier once, the excess proton often immediately recrosses
the barrier. In order to quantify this effect, the transfer paths are grouped
into transfer events. A transfer event is defined by subsequently occurring
transfer paths without recrossing of the same minimum and is sorted by the
number of these crossings. An example for this definition is given in fig.
S12A. The normalized distribution of the transfer events by the number of
recrossings, in the following called the recrossing-number distribution
$p_{\mathrm{RN}}(n)$, is given in fig. S12B. For low barriers up to 40
subsequent recrossings are observed. For higher barriers the distribution is
shifted to lower numbers of recrossings. Nevertheless, for all barrier heights
a significant fraction of transfer events consist of multiple recrossings,
$\sum_{i=1}^{\infty}p_{\mathrm{RN}}(n)=0.47$ for
$R_{\mathrm{OO}}=$2.56\text{\,}\mathrm{\mbox{\AA}}$$, 0.32 for
$R_{\mathrm{OO}}=$2.64\text{\,}\mathrm{\mbox{\AA}}$$ and 0.14 for
$R_{\mathrm{OO}}=$2.72\text{\,}\mathrm{\mbox{\AA}}$$. The mean number of
recrossings is $\sum_{i=0}^{\infty}n\ p_{\mathrm{RN}}(n)=1.84$ for
$R_{\mathrm{OO}}=$2.56\text{\,}\mathrm{\mbox{\AA}}$$, 0.70 for
$R_{\mathrm{OO}}=$2.64\text{\,}\mathrm{\mbox{\AA}}$$ and 0.41 for
$R_{\mathrm{OO}}=$2.72\text{\,}\mathrm{\mbox{\AA}}$$.

\begin{overpic}[width=212.47617pt]{{/../figs/zundel_traj_recrExp}.pdf}
\put(2.0,53.0){\huge\bf A} \put(27.0,58.0){$\overbrace{\quad}^{n=2}$}
\put(37.0,58.0){$\overbrace{\ }^{n=1}$}
\put(62.0,58.0){$\overbrace{\quad\quad\quad\quad\quad}^{n=0}$}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs3/tp/zundel_tp_recrStat}.pdf}
\put(2.0,53.0){\huge\bf B} \end{overpic}

Figure S12: A Example trajectory of the excess-proton coordinate $d(t)$ in the
H5O${}_{2}{}^{+}$ cation with constrained
$R_{\mathrm{OO}}$=$2.64\text{\,}\mathrm{\mbox{\AA}}$, showing nine transfer
paths between the minima of the free energy denoted as grey dotted lines. The
first three transfer paths belong to a single transfer event (with two
recrossings). The subsequent two transfer paths also belong to a single
transfer event (with one recrossing). The remaining four transfer paths are
all defined to not be recrossing. It follows $p_{\mathrm{RN}}(0)=4/6$,
$p_{\mathrm{RN}}(1)=1/6$ and $p_{\mathrm{RN}}(2)=1/6$. B Recrossing-number
probability distributions $p_{\mathrm{RN}}(n)$ normalized such that
$\sum_{n=0}^{\infty}p_{\mathrm{RN}}(n)=1$, for different $R_{\mathrm{OO}}$.

##  VII IR power spectra from linear-response theory

Assuming linear response of an observable $x(t)$ with respect to a force that
couples to an observable $y(t)$, the response function $\chi_{xy}(t)$ is
related to the correlation function $C_{xy}(t^{\prime})=\langle
x(t+t^{\prime})y(t)\rangle$

$\displaystyle\chi_{xy}(t)=-\frac{1}{k_{B}T}\Theta(t)\frac{d}{dt}C_{xy}(t),$
(S9)

where $k_{B}T$ is the thermal energy. Realizing that $\chi(t)$ is single-
sided, the Fourier transform is calculated as

$\displaystyle\widetilde{\chi}_{xy}(\omega)$
$\displaystyle=-\frac{1}{k_{B}T}\int_{0}^{\infty}dt\ e^{i\omega
t}\frac{d}{dt}C_{xy}(t)$
$\displaystyle=-\frac{1}{k_{B}T}\left(C_{xy}(0)-i\omega\int_{0}^{\infty}dt\
e^{i\omega t}C_{xy}(t)\right)$
$\displaystyle=-\frac{1}{k_{B}T}\left(C_{xy}(0)-i\omega\widetilde{C}_{xy}^{+}(\omega)\right),$
(S10)

where the superscript + denotes a single-sided Fourier transform. In case of
$x=y$, $C_{xx}(t)$ is an autocorrelation function, which is real and
symmetric, therefore it follows for the imaginary part of the response
function in Fourier space

$\displaystyle\widetilde{\chi}_{xx}^{\prime\prime}(\omega)$
$\displaystyle=\frac{1}{k_{B}T}\omega\operatorname{Re}(\widetilde{C}^{+}_{xx}(\omega))$
(S11)
$\displaystyle=\frac{1}{k_{B}T}\frac{\omega}{2}\widetilde{C}_{xx}(\omega).$
(S12)

When computing the power spectra of a stochastic process $x(t)$, limited to
the time domain $[0,L_{t}]$, the Wiener-Kintchine theoreme, eq. (S64) in
section XVI, can be used to express $\widetilde{C}_{xx}(\omega)$ in terms of
$\tilde{x}(\omega)$, turning eq. (S12) into

$\displaystyle\widetilde{\chi}_{xx}^{\prime\prime}(\omega)=\frac{\omega}{2k_{B}TL_{t}}|\tilde{x}(\omega)|^{2}.$
(S13)

In case of $x(t)$ being the polarization $\bm{p}(t)$ of the system, which is
coupled to an external electric field $\bm{E}(t)$, the dimensionless
dielectric susceptibility $\chi(t)$ is given by

$\displaystyle\widetilde{\chi}(\omega)=\frac{1}{V\epsilon_{0}l}\langle\widetilde{\chi}_{\bm{p}\bm{p}}(\omega)\rangle,$
(S14)

where $\epsilon_{0}$ is the vacuum permittivity, $V$ is the system volume and
an average is performed over the $l$ dimensions of $\bm{p}$.

##  VIII Spectral cross contributions of excess-proton dynamics in the
H5O${}_{2}{}^{+}$ cation

A decomposition of a trajectory $x(t)$ into two parts $x(t)=x_{1}(t)+x_{2}(t)$
gives rise to three contributions in the total power spectrum

$\displaystyle\omega\widetilde{\chi}_{xx}^{\prime\prime}(\omega)$
$\displaystyle=\frac{\omega^{2}}{2k_{B}T}\left[\widetilde{C}_{1}(\omega)+\widetilde{C}_{2}(\omega)+2\widetilde{C}_{1,2}(\omega)\right]$
$\displaystyle=\omega\left[\widetilde{\chi}^{\prime\prime}_{1}(\omega)+\widetilde{\chi}^{\prime\prime}_{2}(\omega)+\widetilde{\chi}^{\prime\prime}_{1,2}(\omega)\right],$
(S15)

where the cross-correlation contribution $\chi^{\prime\prime}_{1,2}(\omega)$
is defined such that it equals the difference spectrum

$\displaystyle\widetilde{\chi}^{\prime\prime}_{\mathrm{diff}}(\omega)$
$\displaystyle=\widetilde{\chi}_{xx}^{\prime\prime}-\widetilde{\chi}^{\prime\prime}_{1}(\omega)-\widetilde{\chi}^{\prime\prime}_{2}(\omega)=\widetilde{\chi}^{\prime\prime}_{1,2}(\omega)$
$\displaystyle=\frac{\omega}{k_{B}T}\widetilde{C}_{1,2}(\omega).$ (S16)

A positive cross-correlation spectrum hints to in-phase motion, a negative
cross-correlation spectrum to out-of-phase motion of $x_{1}(t+t^{\prime})$ and
$x_{2}(t)$ at a given frequency.

\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tp_decompCrossD_d2.56_log}.pdf}
\put(-2.0,88.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tp_decompCrossD_d2.64_log}.pdf}
\put(-2.0,88.0){\huge\bf B}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tp_decompCrossD_d2.72_log}.pdf}
\put(-2.0,88.0){\huge\bf C} \end{overpic}

Figure S13: Spectral decomposition of the excess-proton dynamics in the
H5O${}_{2}{}^{+}$ cation for different constrained $R_{\mathrm{OO}}$ (A:
$R_{\mathrm{OO}}=$2.56\text{\,}\mathrm{\mbox{\AA}}$$, B:
$R_{\mathrm{OO}}=$2.64\text{\,}\mathrm{\mbox{\AA}}$$, C:
$R_{\mathrm{OO}}=$2.72\text{\,}\mathrm{\mbox{\AA}}$$), similar to the results
in fig. LABEL:decompB in the main text (colored lines). The different cross-
contributions are shown in red and negative cross-contributions as red broken
lines.

Therefore the decomposition of the excess-proton dynamics in the time domain
according to $d(t)=d_{\mathrm{TW}}(t)+d_{\mathrm{TP}}(t)+d_{\mathrm{NM}}(t)$
into a transfer-waiting $d_{\mathrm{TW}}(t)$, a transfer-path
$d_{\mathrm{TP}}(t)$ and a normal-mode contribution $d_{\mathrm{NM}}(t)$ as
introduced in fig. LABEL:decomp, in the main text, produces spectral cross-
contributions. The spectral contributions
$\omega\widetilde{\chi}^{\prime\prime}_{\mathrm{TW}}$,
$\omega\widetilde{\chi}^{\prime\prime}_{\mathrm{TP}}$ and
$\omega\widetilde{\chi}^{\prime\prime}_{\mathrm{NM}}$ as well as the spectral
cross contributions $\omega\widetilde{\chi}^{\prime\prime}_{\mathrm{TW,NM}}$,
$\omega\widetilde{\chi}^{\prime\prime}_{\mathrm{TP,TW}}$ and
$\omega\widetilde{\chi}^{\prime\prime}_{\mathrm{NM,TP}}$ resulting from this
decomposition are shown in fig. S13.

##  IX Dynamics of a two-state process

Consider an arbitrary two-state process $d(t)$ characterized by a a jump
distance $D$, a jump-time probability distribution, which is the transfer-
waiting-time probability distribution $p_{\mathrm{TW}}(t)$ with
$\int_{0}^{\infty}p_{\mathrm{TW}}(t)dt=1$, and a survival probability
$q_{\mathrm{TW}}(t)=\int_{t}^{\infty}p_{\mathrm{TW}}(t^{\prime})dt^{\prime}$.
The conditional expectation values $\langle d(t)\rangle|_{d(0)}$ for the
process starting at time $t=0$ at either $d(0)=0$ or $d(0)=D$ are thus given
as

$\displaystyle\langle d(t)\rangle|_{d(0)=0}$ $\displaystyle=\cancel{0\
p_{00}(t)}+D\ p_{0D}(t),$ (S17) $\displaystyle\langle d(t)\rangle|_{d(0)=D}$
$\displaystyle=\cancel{0\ p_{D0}(t)}+D\ p_{\mathrm{DD}}(t),$ (S18)

where $p_{xy}(t)$ denotes the probability to be at $d(t)=y$ when starting at
$d(0)=x$. In the long-time limit each state has equal probability. The
autocorrelation function, $C(t)=\langle d(t)d(0)\rangle$, is therefore given
as

$\displaystyle\langle d(t)d(0)\rangle$
$\displaystyle=\frac{1}{2}\left(\cancel{0\ \langle d(t)\rangle|_{d(0)=0}}+D\
\langle d(t)\rangle|_{d(0)=D}\right)$ (S19)
$\displaystyle=\frac{D^{2}}{2}p_{\mathrm{DD}}(t).$ (S20)

In order to express the probability $p_{\mathrm{DD}}(t)$ in terms of the
first-passage-time distribution $p_{\mathrm{TW}}(t)$ and survival probability
$q_{\mathrm{TW}}(t)$, all possible jumps within time $t$ have to be
considered, as illustrated in fig. S14

$\displaystyle\begin{split}p_{\mathrm{DD}}(t)=&\int_{0}^{\infty}dt_{D}q_{\mathrm{TW}}(t_{D})\sum_{N=0}^{\infty}\prod_{j=1}^{N}\int_{0}^{\infty}dt_{j}^{D}p_{\mathrm{TW}}(t_{j}^{D})\\\
&\int_{0}^{\infty}dt_{j}^{0}p_{\mathrm{TW}}(t_{j}^{0})\
\delta\left[t-t_{D}-\sum_{j=1}^{N}(t_{j}^{D}+t_{j}^{0})\right].\end{split}$
(S21)

\begin{overpic}[width=151.76964pt]{{/../figs/barrier_theory4}.pdf}
\end{overpic} Figure S14: Schematic of a binary jump process $d(t)$ with jumps
separated by residence times $t_{j}^{0}$ and $t_{j}^{D}$ following the
notation considered in eq.(S21). A doublewell potential, U($d$), is shown to
highlight the relation of the binary jump process to the barrier-crossing
dynamics.

Since $p_{\mathrm{TW}}(t)$ and $q_{\mathrm{TW}}(t)$ are single-sided, a
single-sided Fourier transform is performed,

$\displaystyle\begin{split}\tilde{p}_{\mathrm{DD}}^{+}(\omega)=&\int_{0}^{\infty}dt\
e^{i\omega
t}\int_{0}^{\infty}dt_{D}q_{\mathrm{TW}}(t_{D})\sum_{N=0}^{\infty}\prod_{j=1}^{N}\int_{0}^{\infty}dt_{j}^{D}p_{\mathrm{TW}}(t_{j}^{D})\\\
&\int_{0}^{\infty}dt_{j}^{0}p_{\mathrm{TW}}(t_{j}^{0})\
\delta\left[t-t_{D}-\sum_{j=1}^{N}(t_{j}^{D}+t_{j}^{0})\right]\end{split}$
(S22)
$\displaystyle\begin{split}=&\int_{0}^{\infty}dt_{D}q_{\mathrm{TW}}(t_{D})\sum_{N=0}^{\infty}\prod_{j=1}^{N}\int_{0}^{\infty}dt_{j}^{D}p_{\mathrm{TW}}(t_{j}^{D})\\\
&\int_{0}^{\infty}dt_{j}^{0}p_{\mathrm{TW}}(t_{j}^{0})\
e^{i\omega\left[t_{D}+\sum_{j=1}^{N}(t_{j}^{D}+t_{j}^{0})\right]}\end{split}$
(S23) $\displaystyle=$
$\displaystyle\tilde{q}_{\mathrm{TW}}(\omega)\sum_{N=0}^{\infty}\tilde{p}_{\mathrm{TW}}(\omega)^{2N}$
(S24) $\displaystyle=$
$\displaystyle\frac{\tilde{q}_{\mathrm{TW}}(\omega)}{1-\tilde{p}_{\mathrm{TW}}(\omega)^{2}},$
(S25)

and it follows for the single-sided Fourier transformed autocorrelation

$\displaystyle\widetilde{C}^{+}(\omega)=\frac{D^{2}}{2}\frac{\tilde{q}_{\mathrm{TW}}(\omega)}{1-\tilde{p}_{\mathrm{TW}}(\omega)^{2}}.$
(S26)

From eqs. (S10) and (S26) we obtain our expression for the response function
of the binary jump process

$\displaystyle\widetilde{\chi}_{DD}(\omega)=-\frac{D^{2}}{2k_{B}T}\left(1-\frac{i\omega\tilde{q}_{\mathrm{TW}}(\omega)}{1-\tilde{p}_{\mathrm{TW}}(\omega)^{2}}\right).$
(S27)

In case of the jumping variable $d(t)$ being a polarization with polarization
jump $2qd^{*}_{\rm TW}$, where $q$ is the charge, the dielectric
susceptibility of the binary polarization jump process reads

$\displaystyle\widetilde{\chi}(\omega)=-\frac{2q^{2}{d^{*}_{\rm
TW}}^{2}}{V\epsilon_{0}k_{B}T}\left(1-\frac{i\omega\tilde{q}_{\mathrm{TW}}(\omega)}{1-\tilde{p}_{\mathrm{TW}}(\omega)^{2}}\right),$
(S28)

and eventually its power spectrum, proportional to the imaginary part of the
dielectric susceptibility, is obtained as

$\displaystyle\omega\widetilde{\chi}^{\prime\prime}(\omega)=\frac{2q^{2}{d^{*}_{\rm
TW}}^{2}}{V\epsilon_{0}k_{B}T}\
\text{Re}\left(\frac{\omega^{2}\tilde{q}_{\mathrm{TW}}(\omega)}{1-\tilde{p}_{\mathrm{TW}}(\omega)^{2}}\right).$
(S29)

##  X Power spectrum of the damped harmonic oscillator

The absorbed power $\omega\widetilde{\chi}_{xx}^{\prime\prime}(\omega)$ of the
damped harmonic oscillator described by the differential equation

$\displaystyle m\ddot{x}(t)=-\gamma\dot{x}(t)-kx(t)+F_{\mathrm{ext}}(t),$
(S30)

is computed from the linear response in Fourier space

$\displaystyle\widetilde{\chi}_{xx}(\omega)$
$\displaystyle=\frac{\tilde{x}(\omega)}{\tilde{F}_{\mathrm{ext}}(\omega)}$
(S31) $\displaystyle=(k-m\omega^{2}-i\gamma\omega)^{-1}$ (S32)
$\displaystyle=\frac{k-m\omega^{2}+i\gamma\omega}{(k-m\omega^{2})^{2}+\gamma^{2}\omega^{2}},$
(S33)

where $\tilde{x}(\omega)$ is the oscillating variable, $m$ the mass, $\gamma$
the friction coefficient, $k$ the spring constant of the harmonic potential
and $\tilde{F}_{\mathrm{ext}}(\omega)$ an external force. For the power
spectrum follows

$\displaystyle\omega\widetilde{\chi}_{xx}^{\prime\prime}(\omega)=\frac{\gamma\omega^{2}}{(k-m\omega^{2})^{2}+\gamma^{2}\omega^{2}},$
(S34)

which by introducing the time scales $\tau=2\gamma/k$, $\tau_{m}=\sqrt{m/k}$
and length scale $D$ with $D^{2}=k_{B}T/k$ converts to

$\displaystyle\omega\widetilde{\chi}_{xx}^{\prime\prime}(\omega)=\frac{2D^{2}}{k_{B}T}\frac{\tau\omega^{2}}{4(1-\tau_{m}^{2}\omega^{2})^{2}+\tau^{2}\omega^{2}}.$
(S35)

In case of the oscillating variable $x(t)$ being a polarization with
polarization jump $qD$, where $q$ is the charge, the dielectric susceptibility
reads

$\displaystyle\omega\widetilde{\chi}^{\prime\prime}(\omega)=\frac{2q^{2}D^{2}}{V\epsilon_{0}k_{B}T}\frac{\tau\omega^{2}}{4(1-\tau_{m}^{2}\omega^{2})^{2}+\tau^{2}\omega^{2}}.$
(S36)

In spectroscopy this is known as a Lorentz band shape, which in the overdamped
case, $\tau_{m}\to 0$, reads

$\displaystyle\omega\widetilde{\chi}^{\prime\prime}(\omega)=\frac{2q^{2}D^{2}}{V\epsilon_{0}k_{B}T}\frac{\tau\omega^{2}}{4+\tau^{2}\omega^{2}},$
(S37)

known as the Debye band shape [20].

##  XI Transfer-path shape

Kim and Netz [21] derived the transfer-path-time shape $t(d_{\text{TP}})$ over
a harmonic barrier by an exact calculation, valid for arbitrary barrier
height, as

$\displaystyle t(d_{\text{TP}})$ $\displaystyle=$
$\displaystyle{\tau}-\frac{2\gamma{d^{*}_{\mathrm{TW}}}^{2}}{k_{B}TU}\int_{\sqrt{U}}^{\sqrt{U}(d_{\text{TP}}/d^{*}_{\mathrm{TW}}-1)}\text{d}y\left(\frac{\text{erf}(y)-\text{erf}(\sqrt{U})}{\text{erf}(\sqrt{U}(d_{\text{TP}}/d^{*}_{\mathrm{TW}}-1))-\text{erf}(\sqrt{U})}-\frac{1}{2}\right)D_{+}(y),$
(S38)

where $D_{+}(x)=e^{-x^{2}}\int_{0}^{x}dt\ e^{t^{2}}$ is the Dawson integral
function, $\gamma$ is the friction constant, $U=U_{0}/k_{B}T$ is the
dimensionless barrier height, $d^{*}_{\mathrm{TW}}$ a length scale and $\tau$
a time scale. Note that the shape function is given by time $t$ as a function
of position $d$. The second term in eq. (S38) vanishes for
$d_{\text{TP}}=2d^{*}_{\mathrm{TW}}$ and reduces to $-{\tau}$ for
$d_{\text{TP}}=0$ and therefore $t(d_{\text{TP}}=0)=0$ and
$t(d_{\text{TP}}=2d^{*}_{\mathrm{TW}})={\tau}$. For a variable $s=y/\sqrt{U}$,
eq. (S38) is rewritten as

$\displaystyle t(d_{\text{TP}})$ $\displaystyle=$
$\displaystyle{\tau}-\frac{2\gamma{d^{*}_{\mathrm{TW}}}^{2}}{k_{B}T\sqrt{U}}\int_{1}^{d_{\text{TP}}/d^{*}_{\mathrm{TW}}-1}\text{ds}\left(\frac{\text{erf}(\sqrt{U}s)-\text{erf}(\sqrt{U})}{\text{erf}(\sqrt{U}(d_{\text{TP}}/d^{*}_{\mathrm{TW}}-1))-\text{erf}(\sqrt{U})}-\frac{1}{2}\right)D_{+}(\sqrt{U}s),$
(S39)

For large $U\gg 1$ and
$d^{*}_{\mathrm{TW}}<d_{\text{TP}}<2d^{*}_{\mathrm{TW}}$, eq. (S39) has an
asymptotic expression

$\displaystyle t(d_{\text{TP}})$
$\displaystyle\approx{\tau}+\frac{\gamma{d^{*}_{\mathrm{TW}}}^{2}}{k_{B}T}\int_{1}^{d_{\text{TP}}/d^{*}_{\mathrm{TW}}-1}\text{d}t\frac{D_{+}(\sqrt{U}t)}{\sqrt{U}}$
$\displaystyle={\tau}+\frac{\gamma{d^{*}_{\mathrm{TW}}}^{2}}{2k_{B}TU}\ln(d_{\text{TP}}/d^{*}_{\mathrm{TW}}-1),$
(S40)

where we use $D_{+}(\sqrt{U}t)/\sqrt{U}\approx 1/(2Us)$.

Therefore, further using the symmetric nature of $t(d_{\text{TP}})$ in the
limit $U\rightarrow\infty$, we obtain the asymptotic expression for
$t(d_{\text{TP}})$ as

$\displaystyle
t(d_{\text{TP}})=\begin{cases}-\frac{\gamma{d^{*}_{\mathrm{TW}}}^{2}}{2k_{B}TU}\ln(1-d_{\text{TP}}/d^{*}_{\mathrm{TW}}),\leavevmode\nobreak\
&\text{for}\leavevmode\nobreak\ 0<d_{\text{TP}}<d^{*}_{\mathrm{TW}}\\\
{\tau}+\frac{\gamma{d^{*}_{\mathrm{TW}}}^{2}}{2k_{B}TU}\ln(d_{\text{TP}}/d^{*}_{\mathrm{TW}}-1),\leavevmode\nobreak\
&\text{for}\leavevmode\nobreak\
d^{*}_{\mathrm{TW}}<d_{\text{TP}}<2d^{*}_{\mathrm{TW}}.\end{cases}$ (S41)

It is straightforward to solve eq. (S41) for $d_{\text{TP}}$, yielding

$\displaystyle
d_{\text{TP}}(t)=\begin{cases}-d^{*}_{\mathrm{TW}}\left(e^{-\frac{2k_{B}TUt}{\gamma{d^{*}_{\mathrm{TW}}}^{2}}}-1\right),\leavevmode\nobreak\
&\text{for}\leavevmode\nobreak\ 0<t<{\tau}/2\\\
d^{*}_{\mathrm{TW}}\left(e^{\frac{2k_{B}TU(t-\tau)}{\gamma{d^{*}_{\mathrm{TW}}}^{2}}}+1\right),\leavevmode\nobreak\
&\text{for}\leavevmode\nobreak\ {\tau}/2<t<{\tau}.\end{cases}$ (S42)

Using the curvature parameter
$\kappa=\gamma{d^{*}_{\mathrm{TW}}}^{2}/(2k_{B}TU)$, and shifting the
variables $t\rightarrow t+{\tau}/2$ and $d_{\text{TP}}\rightarrow
d_{\text{TP}}+d^{*}_{\mathrm{TW}}$ to fulfill $d_{\text{TP}}(t=0)=0$, we
arrive at the leading order expression as the sum of the above two

$\displaystyle d_{\text{TP}}(t)$
$\displaystyle=d^{*}_{\mathrm{TW}}\left[e^{-\frac{\tau-2t}{2\kappa}}-e^{-\frac{\tau+2t}{2\kappa}}\right]$
(S43)
$\displaystyle=d^{*}_{\mathrm{TW}}\left[e^{t/\kappa}-e^{-t/\kappa}\right]/e^{\frac{\tau}{2\kappa}}.$

This expression is easily compared to the transfer-path shape (in the presence
of a harmonic potential) derived from the path-integral approach (equivalent
to eq. 3 in the main text) [22, 23]

$\displaystyle
d_{\text{TP}}(t)={d^{*}_{\mathrm{TW}}}\left[e^{t/\kappa}-e^{-t/\kappa}\right]/{\cal
N},$ (S44)

with a slightly different normalization factor ${\cal
N}=e^{\frac{\tau}{2\kappa}}-e^{-\frac{\tau}{2\kappa}}$. Note that the
difference vanishes in the high-barrier limit $U\rightarrow\infty$, i.e.
$\kappa\rightarrow 0$, in which limit eq. (S43) derived from reference [21]
becomes equivalent to eq. (S44) derived from the path-integral approach.

##  XII Spectral signatures of transfer paths

The IR spectral signature of transfer paths is derived by modeling the
Fourier-transformed transfer-path contribution
$\tilde{d}_{\mathrm{TP}}(\omega)$ based on the transfer-path time
$\tau_{\mathrm{TP}}$ and the recrossing-number probability distribution
$p_{\mathrm{RN}}(n)$.

\begin{overpic}[width=212.47617pt]{{/../figs/tp_theory7}.pdf}
\put(-2.0,65.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs3/tp/zundel_tp_d2.56_MinMaxTimes_ExpN}.eps}
\put(-2.0,65.0){\huge\bf B}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs3/tp/zundel_tp_d2.64_MinMaxTimes_ExpN}.eps}
\put(-2.0,65.0){\huge\bf C}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs3/tp/zundel_tp_d2.72_MinMaxTimes_ExpN}.eps}
\put(-2.0,65.0){\huge\bf D}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs3/zundel_delta_fit_tps_short}.pdf}
\put(-2.0,57.0){\huge\bf E}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs3/tp/zundel_tp_recrTimes_ExpN}.pdf}
\put(-2.0,57.0){\huge\bf F} \end{overpic}

Figure S15: A Model for the transfer-path trajectory $d_{\mathrm{TP}}(t)$ with
length $\tau_{\mathrm{TW}}$. In this example the initial transfer path
$d^{0}_{\mathrm{TP}}(t)$ with duration $\tau_{\mathrm{TP}}$ is centered around
$t=0$ and followed by $N=2$ recrossing transfer paths, where $N$ is drawn from
the recrossing-number probability distribution $p_{\mathrm{RN}}(N)$. The
recrossing transfer path are alternating and shifted relative to each other by
the recrossing time $\tau_{\mathrm{RT}}\approx\tau_{\mathrm{TP}}$. B–D
Transfer-path-time probability distributions measured between the turning
points ($p_{\mathrm{TP}}(t)$, B-D), between the minima of the free-energy
landscape, ($p^{\mathrm{w}}_{\mathrm{TP}}(t)$, E) and recrossing-time
distributions measured between recrossing of $d=0$ ($p_{\mathrm{RT}}(t)$, F).
The distributions (data points) are fitted to
$p_{\mathrm{TP}}(t)=\frac{t^{\beta-1}}{(\beta-1)!}\left(\frac{\beta}{\tau}\right)^{\beta}e^{-\beta
t/\tau}$ (solid, broken or dotted lines). Note that only the red curves in B–D
are considered in the calculation of the TP spectral signature.

Recall that the power spectrum $\omega\widetilde{\chi}^{\prime\prime}(\omega)$
of a stochastic process $x(t)$ limited to the time domain $[0,L_{t}]$ is
computed from the Fourier-transformed expressions $\tilde{x}(\omega)$ as
presented in section VII. In the following an expression for
$\tilde{d}_{\mathrm{TP}}(\omega)$ is derived, based on the model for
$d_{\mathrm{TP}}(t)$ illustrated in fig. S15A. The mean transfer path is
expected to repeat on average with period $\tau_{\mathrm{TW}}$. The model is
therefore constrained to $[0,\tau_{TW}]$ and using eq. (S13) the IR power
spectrum of the transfer-path contribution can then be written as

$\displaystyle\omega\widetilde{\chi}_{\mathrm{TP}}^{\prime\prime}(\omega)=\omega^{2}\frac{q^{2}}{V\epsilon_{0}k_{B}T\tau_{\mathrm{TW}}}\tilde{d}_{\mathrm{TP}}(\omega)\tilde{d}_{\mathrm{TP}}^{*}(\omega).$
(S45)

$d_{\mathrm{TP}}(t)$ is modeled by a single mean transfer path with shape
$d^{0}_{\mathrm{TP}}(t)$ that is followed by a subsequent number of recrossing
transfer paths, according to the recrossing-number probability distribution
$p_{\mathrm{RN}}(n)$

$\displaystyle d_{\mathrm{TP}}(t)$
$\displaystyle=\sum_{n=0}^{\infty}p_{\mathrm{RN}}(n)\sum_{m=0}^{n}\
(-1)^{m}d^{0}_{\mathrm{TP}}\left(t-m\tau_{\text{TP}}\right),$ (S46)

where the factor $(-1)^{m}$ accounts for the alternation of recrossing
transfer paths that are going up and down. A Fourier transform with respect to
$t$ turns the previous eq. (S46) into

$\displaystyle\tilde{d}_{\mathrm{TP}}(\omega)$
$\displaystyle=\sum_{n=0}^{\infty}p_{\mathrm{RN}}(n)\sum_{m=0}^{n}\
(-1)^{m}\int_{-\infty}^{\infty}dt\ e^{i\omega t}\
d^{0}_{\mathrm{TP}}(t-m\tau_{\text{TP}})$
$\displaystyle=\tilde{d}^{0}_{\mathrm{TP}}(\omega)\
\sum_{n=0}^{\infty}p_{\mathrm{RN}}(n)\sum_{m=0}^{n}\ (-1)^{m}e^{-i\omega
m\tau_{\text{TP}}}.$ (S47)

\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tpFitMinMaxAvgD_d2.56}.png}
\put(-2.0,55.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tpFitMinMaxAvgD_d2.64}.png}
\put(-2.0,55.0){\huge\bf B}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tpFitMinMaxAvgD_d2.72}.png}
\put(-2.0,55.0){\huge\bf C} \end{overpic}

Figure S16: Ensembles of transfer paths sampled from the trajectories of the
excess proton in the H5O${}_{2}{}^{+}$ cation for different constrained oxygen
distances $R_{\mathrm{OO}}$ (A:
$R_{\mathrm{OO}}=$2.56\text{\,}\mathrm{\mbox{\AA}}$$, B:
$R_{\mathrm{OO}}=$2.64\text{\,}\mathrm{\mbox{\AA}}$$, C:
$R_{\mathrm{OO}}=$2.72\text{\,}\mathrm{\mbox{\AA}}$$), scaled to their
individual transfer-path times $t_{\mathrm{TP}}$ between the respective
turning points. The blue dashed lines indicate the mean initial and final
values, which are used to estimate the parameter $d^{*}_{\mathrm{TP}}$, the
mean transfer-path distance along $d$. The horizontal black dashed lines
indicate the minima of the free energy, which are used to estimate the
parameter $d^{*}_{\mathrm{TW}}$. Mean transfer paths between the minima of the
free energy (colored lines) are fitted using eq. (LABEL:eq:tpPathTheo). The
fits are shown as red dotted lines. Fit parameters are given in the legend.

An expression for the mean transfer-path shape $d^{0}_{\mathrm{TP}}(t)$ is
obtained by regarding the ensembles of rescaled transfer paths, which are
shown in fig. S16 for various $R_{\mathrm{OO}}$, along with the mean transfer
paths, obtained from space-averaging the ensembles of rescaled transfer paths
at each rescaled time between between the respective turning points of the
trajectory. The mean value at the turning points defines the length scale
$d^{*}_{\rm TP}$. The single mean transfer path $d^{0}_{\mathrm{TP}}(t)$
reaching between the turning points is then modeled by a truncated straight
line

$\displaystyle d^{0}_{\mathrm{TP}}(t)$ $\displaystyle=2d^{*}_{\rm
TP}\frac{t}{\tau_{\text{TP}}}\left(\theta\left(t\right)-\theta\left(t-\tau_{\text{TP}}\right)\right),$
(S48)

with the Fourier transform

$\displaystyle\tilde{d}^{0}_{\mathrm{TP}}(\omega)$ $\displaystyle=d^{*}_{\rm
TP}\frac{e^{i\omega\tau_{\text{TP}}}\left(2-i\omega\tau_{\text{TP}}\right)-2-i\omega\tau_{\text{TP}}}{\omega^{2}\tau_{\text{TP}}}.$
(S49)

The time scale $\tau_{\text{TP}}$ in eq. (S47) is estimated directly from the
simulation data, using the distributions, shown in fig. S15B – D. The
distributions of transfer-path times considering only multiple transfers in
fig. S15B–D are sharply peaked at roughly the same value for all systems
$\tau_{\mathrm{TP}}=$0.023\text{\,}\mathrm{p}\mathrm{s}$$.

To account for subsequent recrossings, the expression eq. (S47) is evaluated
for an exponential decay of the recrossing-number probability distribution,
$p_{\mathrm{RN}}(n)=(1-e^{-\alpha})e^{-\alpha n}$, with decay parameter
$\alpha$ and shown in fig. S17 for the given systems. Using this fit function
the final expression for the transfer-path spectral contribution evaluates to

$\displaystyle\omega\widetilde{\chi}_{\mathrm{TP}}^{\prime\prime}(\omega)=\frac{{2d^{*}_{\rm
TP}}^{2}q^{2}}{V\epsilon_{0}k_{B}T\tau_{\mathrm{TW}}}\frac{e^{\alpha}\left(\omega\tau_{\text{TP}}\cos\left(\frac{\omega\tau_{\text{TP}}}{2}\right)-2\sin\left(\frac{\omega\tau_{\text{TP}}}{2}\right)\right){}^{2}}{\omega^{2}\tau_{\text{TP}}^{2}\left(\cosh(\alpha)+\cos\left(\omega\tau_{\text{TP}}\right)\right)}.$
(S50)

\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tp_recrStat_expFit_d2.56}.pdf}
\put(-2.0,55.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tp_recrStat_expFit_d2.64}.pdf}
\put(-2.0,55.0){\huge\bf B}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tp_recrStat_expFit_d2.72}.pdf}
\put(-2.0,55.0){\huge\bf C} \end{overpic}

Figure S17: Recrossing-number probability distributions $p_{\mathrm{RN}}(n)$
for different constrained oxygen distances $R_{\mathrm{OO}}$ (A:
$R_{\mathrm{OO}}=$2.56\text{\,}\mathrm{\mbox{\AA}}$$, B:
$R_{\mathrm{OO}}=$2.64\text{\,}\mathrm{\mbox{\AA}}$$, C:
$R_{\mathrm{OO}}=$2.72\text{\,}\mathrm{\mbox{\AA}}$$), normalized to
$\sum_{n=0}^{\infty}p_{\mathrm{RN}}=1$ and fits according to
$p_{\mathrm{RN}}(n)=(1-e^{-\alpha})e^{-\alpha n}$. The grey shaded areas show
the variation of $\alpha\pm 20\%$

An approximate expression for eq. (S50) is derived by factorizing
$\omega\widetilde{\chi}_{\mathrm{TP}}^{\prime\prime}(\omega)$ into the
recrossing contribution $X_{rec}(\omega,\alpha,\tau_{\text{TP}})$ and the
shape contribution
$|\tilde{d}^{0}_{\mathrm{TP}}(\omega,\tau_{\text{TP}})|^{2}$

$\displaystyle|\tilde{d}^{0}_{\mathrm{TP}}(\omega,\tau_{\text{TP}})|^{2}$
$\displaystyle={d^{*}_{\rm
TP}}^{2}\frac{4\left(\omega\tau_{\text{TP}}\cos\left(\frac{\omega\tau_{\text{TP}}}{2}\right)-2\sin\left(\frac{\omega\tau_{\text{TP}}}{2}\right)\right){}^{2}}{\omega^{4}\tau_{\text{TP}}^{2}},$
(S51) $\displaystyle X_{rec}(\omega,\alpha,\tau_{\text{TP}})$
$\displaystyle=\frac{e^{\alpha}}{2\left(\cosh(\alpha)+\cos\left(\omega\tau_{\text{TP}}\right)\right)},$
(S52)
$\displaystyle\omega\widetilde{\chi}_{\mathrm{TP}}^{\prime\prime}(\omega)$
$\displaystyle=\frac{q^{2}\omega^{2}}{V\epsilon_{0}k_{B}T\tau_{\mathrm{TW}}}X_{rec}(\omega,\alpha,\tau_{\text{TP}})|\tilde{d}^{0}_{\mathrm{TP}}(\omega,\tau_{\text{TP}})|^{2}.$
(S53)

The relevant maximum of $X_{rec}(\omega,\alpha,\tau_{\text{TP}})$ resides at
$\omega=\pi/\tau_{\text{TP}}$, which minimizes the denominator in eq.(S52). A
Taylor expansion of the $\cos$ function in eq.(S52) to second order in
$\omega$ around $\pi/\tau_{\text{TP}}$ leads to

$\displaystyle X_{rec}(\omega,\alpha,\tau_{\text{TP}})$
$\displaystyle=\frac{e^{\alpha}}{2\cosh(\alpha)-2+\left(\omega\tau_{\text{TP}}-\pi\right){}^{2}}.$
(S54)

Furthermore the shape contribution
$|\tilde{d}^{0}_{\mathrm{TP}}(\omega,\tau_{\text{TP}})|^{2}$ can be estimated
around $\omega=\pi/\tau_{\text{TP}}$ by the following expression, which is in
good agreement with the local series expansion

$\displaystyle|\tilde{d}^{0}_{\mathrm{TP}}(\omega,\tau_{\text{TP}})|^{2}$
$\displaystyle={d^{*}_{\rm
TP}}^{2}\frac{64\omega^{2}\tau_{\text{TP}}^{4}}{\pi^{4}\left(\omega\tau_{\text{TP}}+\pi\right){}^{2}},$
(S55)
$\displaystyle\omega\widetilde{\chi}_{\mathrm{TP}}^{\prime\prime}(\omega)$
$\displaystyle=\frac{{d^{*}_{\rm
TP}}^{2}q^{2}}{V\epsilon_{0}k_{B}T\tau_{\mathrm{TW}}}\frac{64e^{\alpha}\omega^{4}\tau_{\text{TP}}^{4}}{\pi^{4}\left(\omega\tau_{\text{TP}}+\pi\right){}^{2}\left(2\cosh(\alpha)-2+\left(\omega\tau_{\text{TP}}-\pi\right){}^{2}\right)}.$
(S56)

Note that eq. (S56) is a multiplication of a function of Debye-shape type with
the shoulder frequency at $\omega=\pi/\tau{{}_{\text{TP}}}$, stemming from the
shape contribution eq. (S55) and a function of Lorentz-shape type with the
resonant frequency at $\omega=\pi/\tau{{}_{\text{TP}}}$, stemming from the
recrossing contribution eq. (S54).

\begin{overpic}[width=212.47617pt]{{/../figs3/barrier/zundel_spectra_tpTheo_lin_all}.png}
\put(-2.0,55.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs3/barrier/zundel_spectra_tpTheo_cos_all}.png}
\put(-2.0,55.0){\huge\bf B} \end{overpic}

Figure S18: Models for the transfer-path spectral contribution given by eq.
(S50) for the linear transfer-path shape (A) and eq. (S59) for the cosine
transfer-path shape (B). The approximations are given by eqs. (S56) and (S60)
and are shown as black dashed lines in A and B. The grey shaded areas show the
variation of $\alpha\pm 20\%$.

If the transfer-path shape $d^{0}_{\mathrm{TP}}(t)$ is alternatively modeled
by a truncated cosine wave with period $2\tau_{\text{TP}}$

$\displaystyle d^{0}_{\mathrm{TP,cos}}(t)$ $\displaystyle=d^{*}_{\rm
TP}\cos\left(\frac{\pi
t}{\tau_{\text{TP}}}\right)\left(\theta\left(t\right)-\theta\left(t-\tau_{\text{TP}}\right)\right),$
(S57)

with the Fourier transform

$\displaystyle\tilde{d}^{0}_{\mathrm{TP,cos}}(\omega)$
$\displaystyle=-d^{*}_{\rm
TP}\frac{i\omega\tau_{\text{TP}}^{2}\left(1+e^{i\omega\tau_{\text{TP}}}\right)}{\omega^{2}\tau_{\text{TP}}^{2}-\pi^{2}},$
(S58)

a slightly different result is obtained.

The final expression

$\displaystyle\omega\widetilde{\chi}_{\mathrm{TP,cos}}^{\prime\prime}(\omega)=\frac{{d^{*}_{\rm
TP}}^{2}q^{2}}{V\epsilon_{0}k_{B}T\tau_{\mathrm{TW}}}\frac{e^{\alpha}\omega^{4}\tau_{\text{TP}}^{4}\left(\cos\left(\omega\tau_{\text{TP}}\right)+1\right)}{\left(\pi^{2}-\omega^{2}\tau_{\text{TP}}^{2}\right){}^{2}\left(\cosh(\alpha)+\cos\left(\omega\tau_{\text{TP}}\right)\right)},$
(S59)

is Taylor-approximated as above to give

$\displaystyle\omega\widetilde{\chi}_{\mathrm{TP,cos}}^{\prime\prime}(\omega)$
$\displaystyle=\frac{{d^{*}_{\rm
TP}}^{2}q^{2}}{V\epsilon_{0}k_{B}T\tau_{\mathrm{TW}}}\frac{e^{\alpha}\omega^{4}\tau_{\text{TP}}^{4}}{\left(\omega\tau_{\text{TP}}+\pi\right){}^{2}\left(2\cosh(\alpha)-2+\left(\omega\tau_{\text{TP}}-\pi\right){}^{2}\right)},$
(S60)

which are both shown in fig. S18B. In fact both approximate equations differ
only by a factor of $64/\pi^{4}\approx 2/3$. The linear shape was chosen for
the presentation in the main text because of the good qualitative agreement of
the mean TP shape in fig. S16 and quantitative agreement of the spectral shape
in fig. S18A.

##  XIII IR spectral decomposition of the excess-proton dynamics in the
constrained H5O2+ cation

Spectral decompositions of the excess-proton dynamics in the H5O${}_{2}{}^{+}$
cation for different constrained $R_{\mathrm{OO}}$ are shown in fig. S19.
Theoretical spectra are shown for the barrier-crossing model derived in
section IX (respective eq. (LABEL:BarrierCrossingSpectrum) in the main text)
and the transfer-path model derived in section XII (respective eq.
(LABEL:TransferPathSpectrum) in the main text), using the distance between the
minima of the free energy in fig. LABEL:systemsB in the main text and fits
applied to the distributions in fig. LABEL:systemsC in the main text and in
section XII.

\begin{overpic}[width=212.47617pt]{{/../figs3/zundel_decompAllDtheo_d2.56_log}.pdf}
\put(-2.0,88.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=212.47617pt]{{/../figs3/zundel_decompAllDtheo_d2.72_log}.pdf}
\put(-2.0,88.0){\huge\bf B} \end{overpic}

Figure S19: Spectral decomposition of the excess-proton dynamics in the
H5O${}_{2}{}^{+}$ cation for different constrained $R_{\mathrm{OO}}$ (A:
$R_{\mathrm{OO}}=$2.56\text{\,}\mathrm{\mbox{\AA}}$$, B:
$R_{\mathrm{OO}}=$2.72\text{\,}\mathrm{\mbox{\AA}}$$) in analogy to the
results in fig. LABEL:decompB in the main text (colored lines). Theoretical
spectra for the models discussed in the previous sections are shown as grey
broken lines.

##  XIV IR spectra of the deuterated H5O2+ cation

\begin{overpic}[width=433.62pt]{{/../figs/zundel_compare_spectra_d_xyz_all}.pdf}
\put(2.0,65.0){\huge\bf A} \put(37.0,65.0){\huge\bf B}
\put(67.0,65.0){\huge\bf C} \put(2.0,45.0){\huge\bf D}
\put(37.0,45.0){\huge\bf E} \put(67.0,45.0){\huge\bf F}
\put(2.0,25.0){\huge\bf G} \put(37.0,25.0){\huge\bf H}
\put(67.0,25.0){\huge\bf I} \end{overpic} Figure S20: Comparison of the IR
spectra of the H5O${}_{2}{}^{+}$ cation with an excess proton (1H+) or an
excess deuteron (2H+) in different directions A, D, G: along the $x$-axis
connecting the two oxygens, B, E, H: along the $yz$-plane C, F, I: isotropic
spectrum. Each row corresponds to a system with a distinct constrained oxygen
distance $R_{\mathrm{OO}}$ given in the first legend. The HOH-bending mode and
OH-stretching mode of the water molecules are indicated in D and E.

In order to test isotope effects, simulations were performed with the excess
proton replaced by an excess deuteron. The resulting IR spectra for various
$R_{\mathrm{OO}}$ and directions are shown in fig. S20 and compared to the
excess-proton data. The regime between
$400\text{\,}\mathrm{c}\mathrm{m}^{-1}2000\text{\,}\mathrm{c}\mathrm{m}^{-1}$
shows various shifts to lower frequencies and also splittings. This H/D
isotope effect is expected since this region is linked to the excess-proton
motion. As expected, the OH-stretching vibrations are not effected. These
vibrations are associated with the flanking water molecules, that are not
deuterated. Fig. S21 compares the H/D isotope effect for the decomposed
spectra. In contrast to the normal-modes involving the excess-proton as well
as the transfer-path signature, which are all shifted by deuteration, the low-
frequency transfer-waiting shoulder does not shift.

\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tp_decompD_HD_d2.56_log}.pdf}
\put(-2.0,88.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tp_decompD_HD_d2.64_log}.pdf}
\put(-2.0,88.0){\huge\bf B}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/tp/zundel_tp_decompD_HD_d2.72_log}.pdf}
\put(-2.0,88.0){\huge\bf C} \end{overpic}

Figure S21: Comparison of the IR spectra of the H5O${}_{2}{}^{+}$ cation with
an excess proton (1H+) or an excess deuteron (2H+) for different
$R_{\mathrm{OO}}$ (A: $R_{\mathrm{OO}}=$2.56\text{\,}\mathrm{\mbox{\AA}}$$, B:
$R_{\mathrm{OO}}=$2.64\text{\,}\mathrm{\mbox{\AA}}$$, C:
$R_{\mathrm{OO}}=$2.72\text{\,}\mathrm{\mbox{\AA}}$$), decomposed as described
in the main text.

##  XV Deuteron transfer-waiting rates

\begin{overpic}[width=138.76157pt]{{/../figs3/d/rates/kramers_d2.56}.pdf}
\put(-2.0,65.0){\huge\bf A}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/d/rates/kramers_d2.64}.pdf}
\put(-2.0,65.0){\huge\bf B}
\end{overpic}\begin{overpic}[width=138.76157pt]{{/../figs3/d/rates/kramers_d2.72}.pdf}
\put(-2.0,65.0){\huge\bf C} \end{overpic}

Figure S22: The transfer-waiting times $\tau_{\mathrm{TW}}$ of the excess
proton (1H+) and the excess deuteron (2H+) of the H5O${}_{2}{}^{+}$ cation for
different constrained oxygen distances $R_{\mathrm{OO}}$ (dotted lines, A:
$R_{\mathrm{OO}}=$2.56\text{\,}\mathrm{\mbox{\AA}}$$, B:
$R_{\mathrm{OO}}=$2.64\text{\,}\mathrm{\mbox{\AA}}$$, C:
$R_{\mathrm{OO}}=$2.72\text{\,}\mathrm{\mbox{\AA}}$$) compared to the
heuristic formula, given by eq. (S61), (straight and broken lines) as a
function of the friction constant $\gamma$. Parameters entering the heuristic
formula are given in the legends.

As discussed in the previous section, no isotope effect is found for the low-
frequency transfer-waiting shoulder of the IR spectrum, i.e. the mean
transfer-waiting time $\tau_{\mathrm{TW}}$ does not depend on the mass of the
reaction coordinate, but solely on friction constant $\gamma$. This is
expected to be the same for the system with either an excess proton or an
excess deuteron. The determination of the friction constant $\gamma$ is not
straight-forward. In fig. S22 the mean transfer-waiting times
$\tau_{\mathrm{TW}}$ measured from the simulations are compared to the
heuristic formula [24]

$\displaystyle\tau_{\mathrm{TW}}=e^{\frac{U_{0}}{k_{B}T}}\left[\frac{m}{\gamma}\frac{k_{B}T}{U_{0}}+\frac{\pi}{2\sqrt{2}}\frac{\gamma{d^{*}_{\rm
TW}}^{2}}{U_{0}}+2\sqrt{\frac{m{d^{*}_{\rm TW}}^{2}}{U_{0}}}\right].$ (S61)

Since the mass $m$, barrier height $U_{0}$ and widths of the barrier
$2{d^{*}_{\rm TW}}$ are known (given in fig. S22), estimates for the friction
constant $\gamma$ are obtained by comparing the simulation data to the
heuristic formula. The obtained values,
$\gamma\approx$400\text{\,}\mathrm{u}\mathrm{/}\mathrm{p}\mathrm{s}$$ for
$R_{\mathrm{OO}}=$2.64\text{\,}\mathrm{\mbox{\AA}}$$ and
$\gamma\approx$900\text{\,}\mathrm{u}\mathrm{/}\mathrm{p}\mathrm{s}$$ for
$R_{\mathrm{OO}}=$2.72\text{\,}\mathrm{\mbox{\AA}}$$, vary greatly and also
deviate from the friction constants fitted to the line-broadening of the
normal-modes, $\gamma=$16\text{\,}\mathrm{u}\mathrm{/}\mathrm{p}\mathrm{s}$$.
This discrepancy highlights the complex frequency dependence of dielectric
friction, as previously discussed by Sedlmeier et al. [25], Brünig et al.
[26]. Note that the theoretical predictions are only valid for
$U_{0}>2\,k_{B}T$, and can therefore not be used for interpretation of the
data for the lowest barrier considered here [24]. Nevertheless, fig. S22
validates the observation of the negligible isotope effect on the mean
transfer-waiting time, since the heuristic formula of eq. (S61) shows little
mass dependence in the predicted regime of $\gamma$, as seen in figs. S22B and
S22C from the negligible difference between the curves for
$m=$1\text{\,}\mathrm{u}$$ (solid lines) and $m=$2\text{\,}\mathrm{u}$$
(dashed lines) in the regime of the measured mean transfer-waiting time
(dotted lines).

##  XVI Wiener-Kintchine theorem

The correlation function $C_{xy}(t)$ of two stochastic processes $x(t)$ and
$y(t)$ limited to the interval $[0,L_{t}]$ is efficiently computed from the
Fourier-transformed expressions $\tilde{x}(\omega)$ and $\tilde{y}(\omega)$
according to

$\displaystyle
C_{xy}(t)=\frac{1}{2\pi(L_{t}-t)}\int_{-\infty}^{\infty}d\omega\ e^{-i\omega
t}\tilde{x}(\omega)\tilde{y}^{*}(\omega),$ (S62)

where the asterix denotes the conjugate form. This is known as the Wiener-
Kintchine theoreme [27]. Both sides of eq. (S62) are Fourier-transformed to
give

$\displaystyle\int_{-\infty}^{\infty}dt\ e^{i\omega t}\
L_{t}\left(1-\frac{t}{L_{t}}\right)C_{xy}(t)$
$\displaystyle=\tilde{x}(\omega)\tilde{y}^{*}(\omega),$ (S63)

which in the limit of large $L_{t}$ reduces to

$\displaystyle\widetilde{C}_{xy}(\omega)$
$\displaystyle=L_{t}^{-1}\tilde{x}(\omega)\tilde{y}^{*}(\omega).$ (S64)

Eq. (S62) can be derived starting off with the definition of the correlation
function

$\displaystyle C_{xy}(t)=\frac{1}{L_{t}-t}\int_{0}^{L_{t}-t}dt^{\prime}\
x(t^{\prime}+t)y(t^{\prime}).$ (S65)

Defining $x(t),y(t)=0$ for $t\not\in[0,L_{t}]$, the integral bounds can
formally be extended

$\displaystyle C_{xy}(t)=\frac{1}{L_{t}-t}\int_{-\infty}^{\infty}dt^{\prime}\
x(t^{\prime}+t)y(t^{\prime}),$ (S66)

and making use of the convolution theorem

$\displaystyle C_{xy}(t)$
$\displaystyle=\frac{1}{4\pi^{2}(L_{t}-t)}\int_{-\infty}^{\infty}dt^{\prime}\int_{-\infty}^{\infty}d\omega\
e^{-i\omega(t+t^{\prime})}\tilde{x}(\omega)\int_{-\infty}^{\infty}d\omega^{\prime}\
e^{-i\omega^{\prime}t^{\prime}}\tilde{y}(\omega^{\prime})$
$\displaystyle=\frac{1}{4\pi^{2}(L_{t}-t)}\int_{-\infty}^{\infty}d\omega\
e^{-i\omega t}\tilde{x}(\omega)\int_{-\infty}^{\infty}d\omega^{\prime}\
\tilde{y}(\omega^{\prime})\int_{-\infty}^{\infty}dt^{\prime}e^{-it^{\prime}(\omega+\omega^{\prime})}$
$\displaystyle=\frac{1}{4\pi^{2}(L_{t}-t)}\int_{-\infty}^{\infty}d\omega\
e^{-i\omega t}\tilde{x}(\omega)\int_{-\infty}^{\infty}d\omega^{\prime}\
\tilde{y}(\omega^{\prime})2\pi\delta(\omega+\omega^{\prime})$
$\displaystyle=\frac{1}{2\pi(L_{t}-t)}\int_{-\infty}^{\infty}d\omega\
e^{-i\omega t}\tilde{x}(\omega)\tilde{y}(-\omega),$ (S67)

noting that $\tilde{y}(-\omega)=\tilde{y}^{*}(\omega)$ for a real function
$y(t)$ in order to obtain eq. (S62).

AIMD

    ab initio molecular-dynamics
DFT

    density functional theory
FPT

    first-passage time
GLE

    generalized Langevin equation
GH

    Grote-Hynes
IR

    infrared
LE

    Langevin equation
MD

    molecular dynamics
MFPT

    mean first-passage time
MFP

    mean first-passage
MSD

    mean squared displacement
NM

    normal-mode
PTP

    $p(\text{TP}|q)$
PME

    particle mesh Ewald [pronk2013gromacs]
PMF

    potential of mean force
PGH

    Pollak-Grabert-Hanggi
RC

    reaction coordinate
RDF

    radial distribution function
RTT

    round-trip time
TP

    transfer path

## References

  * Asmis et al. [2003] Asmis, K. R. _et al._ Gas-Phase Infrared Spectrum of the Protonated Water Dimer.  _Science_ 299, 1375 (2003).
  * Fridgen et al. [2004] Fridgen, T. D., McMahon, T. B., MacAleese, L., Lemaire, J. & Maitre, P.  Infrared spectrum of the protonated water dimer in the gas phase.  _J. Phys. Chem. A_ 108, 9008 (2004).
  * Diken et al. [2005] Diken, E. G. _et al._ Fundamental excitations of the shared proton in the H3O${}_{2}{}^{-}$ and H5O${}_{2}{}^{+}$ complexes.  _J. Phys. Chem. A_ 109, 1487 (2005).
  * Hammer et al. [2005] Hammer, N. I. _et al._ The vibrational predissociation spectra of the H5O${}_{2}{}^{+}\cdot$RGn(RG=Ar,Ne) clusters: Correlation of the solvent perturbations in the free OH and shared proton transitions of the Zundel ion.  _J. Chem. Phys._ 122, 244301 (2005).
  * Dahms et al. [2017] Dahms, F., Fingerhut, B. P., Nibbering, E. T., Pines, E. & Elsaesser, T.  Large-amplitude transfer motion of hydrated excess protons mapped by ultrafast 2D IR spectroscopy.  _Science_ 357, 491 (2017).
  * Sauer and Döbler [2005] Sauer, J. & Döbler, J.  Gas-phase infrared spectrum of the protonated water dimer: Molecular dynamics simulation and accuracy of the potential energy surface.  _ChemPhysChem_ 6, 1706 (2005).
  * Vendrell et al. [2007] Vendrell, O., Gatti, F. & Meyer, H. D.  Full dimensional (15-dimensional) quantum-dynamical simulation of the protonated water dimer. II. Infrared spectrum and vibrational dynamics.  _J. Chem. Phys._ 127, 184303 (2007).
  * Park et al. [2007] Park, M., Shin, I., Singh, N. J. & Kim, K. S.  Eigen and Zundel forms of small protonated water clusters: Structures and infrared spectra.  _J. Phys. Chem. A_ 111, 10692 (2007).
  * Kaledin et al. [2009] Kaledin, M., Kaledin, A. L., Bowman, J. M., Ding, J. & Jordan, K. D.  Calculation of the vibrational spectra of H5O${}_{2}{}^{+}$ and its deuterium-substituted isotopologues by molecular dynamics simulations.  _J. Phys. Chem. A_ 113, 7671 (2009).
  * Huang et al. [2005] Huang, X., Braams, B. J. & Bowman, J. M.  Ab initio potential energy and dipole moment surfaces for H5O${}_{2}{}^{+}$.  _J. Chem. Phys._ 122, 044308 (2005).
  * Baer et al. [2010] Baer, M., Marx, D. & Mathias, G.  Theoretical messenger spectroscopy of microsolvated hydronium and Zundel cations.  _Angew. Chemie - Int. Ed._ 49, 7346 (2010).
  * Agostini et al. [2011] Agostini, F., Vuilleumier, R. & Ciccotti, G.  Infrared spectroscopy of small protonated water clusters at room temperature: An effective modes analysis.  _J. Chem. Phys._ 134, 084302 (2011).
  * Kulig and Agmon [2013] Kulig, W. & Agmon, N.  A ’clusters-in-liquid’ method for calculating infrared spectra identifies the proton-transfer mode in acidic aqueous solutions.  _Nat. Chem._ 5, 29 (2013).
  * Rossi et al. [2014] Rossi, M., Ceriotti, M. & Manolopoulos, D. E.  How to remove the spurious resonances from ring polymer molecular dynamics.  _J. Chem. Phys._ 140, 234116 (2014).
  * Biswas et al. [2017] Biswas, R., Carpenter, W., Fournier, J. A., Voth, G. A. & Tokmakoff, A.  IR spectral assignments for the hydrated excess proton in liquid water.  _J. Chem. Phys._ 146, 154507 (2017).
  * Headrick et al. [2004] Headrick, J. M., Bopp, J. C. & Johnson, M. A.  Predissociation spectroscopy of the argon-solvated H5O${}_{2}{}^{+}$ "zundel" cation in the 1000-1900 cm -1 region.  _J. Chem. Phys._ 121, 11523 (2004).
  * Thämer et al. [2015] Thämer, M., De Marco, L., Ramasesha, K., Mandal, A. & Tokmakoff, A.  Ultrafast 2D IR spectroscopy of the excess proton in liquid water.  _Science_ 350, 78 (2015).
  * Fournier et al. [2018] Fournier, J. A., Carpenter, W. B., Lewis, N. H. & Tokmakoff, A.  Broadband 2D IR spectroscopy reveals dominant asymmetric H5O${}_{2}{}^{+}$ proton hydration structures in acid solutions.  _Nat. Chem._ 10, 932 (2018).
  * Wolke et al. [2016] Wolke, C. T. _et al._ Spectroscopic snapshots of the proton-transfer mechanism in water.  _Science_ 354, 1131 (2016).
  * Debye [1929] Debye, P., _Polar Molecules_ (Chemical Catalog Company, Incorporated, 1929).
  * Kim and Netz [2015] Kim, W. K. & Netz, R. R.  The mean shape of transition and first-passage paths.  _J. Chem. Phys._ 143, 224108 (2015).
  * Faccioli et al. [2006] Faccioli, P., Sega, M., Pederiva, F. & Orland, H.  Dominant pathways in protein folding.  _Phys. Rev. Lett._ 97, 108101 (2006).
  * Cossio et al. [2018] Cossio, P., Hummer, G. & Szabo, A.  Transition paths in single-molecule force spectroscopy.  _J. Chem. Phys._ 148, 123309 (2018).
  * Kappler et al. [2019] Kappler, J., Hinrichsen, V. B. & Netz, R. R.  Non-Markovian barrier crossing with two-time-scale memory is dominated by the faster memory component.  _Eur. Phys. J. E_ 42, 119 (2019).
  * Sedlmeier et al. [2014] Sedlmeier, F., Shadkhoo, S., Bruinsma, R. & Netz, R. R.  Charge/mass dynamic structure factors of water and applications to dielectric friction and electroacoustic conversion.  _J. Chem. Phys._ 140, 054512 (2014).
  * Brünig et al. [2022] Brünig, F. N., Geburtig, O., von Canal, A., Kappler, J. & Netz, R. R.  Time-Dependent Friction Effects on Vibrational Infrared Frequencies and Line Shapes of Liquid Water.  _J. Phys. Chem. B_ 126, 1579 (2022).
  * Wiener [1930] Wiener, N.  Generalized harmonic analysis.  _Acta Math._ 55, 117 (1930).

  *[AIMD]: ab initio molecular-dynamics
  *[IR]: infrared
  *[MD]: molecular dynamics