Decoding the 2D IR spectrum of the aqueous proton with high-level VSCF/VCI calculations

The aqueous proton is one of the cornerstone entities in aqueous chemistry. Beyond standard acid-base chemistry,^1,2 the aqueous proton also plays a key role in catalysis,^3–5 biological energy conversion,^6,7 and hydrogen fuel cell technology.^8,9 Despite its importance, a clear molecular description of the aqueous proton has remained elusive for over a century due to the complexity of the proton’s hydration structures and transport in solution.^2,10–15 Ultrafast fluctuations of the hydrogen bond (H-bond) network in liquid water^16–18 and the low barriers separating different proton hydration configurations are responsible for the complexity of the problem.^19–21 Experimental determination of the distribution of proton hydration configurations would be invaluable as a foundation for studying proton transport in experiments and simulations.

The aqueous proton has been described in terms of two limiting structural motifs, the Eigen¹¹ and Zundel¹² cations. In the Eigen cation H₉O₄⁺, a central hydronium ion is tightly H-bonded to three equivalent waters, whereas the Zundel cation H₅O₂⁺ consists of a central excess proton equally solvated by two flanking waters. These two cations are highly symmetric in the gas phase,^22,23 but the wide variety of aqueous H-bonding configurations in the liquid phase suggest that there is a distribution of H⁺(aq) hydration configurations,^24,25 breaking the symmetries inherent to the Eigen and Zundel complexes. While some simulations have described an aqueous proton structure as a distribution with Eigen and Zundel motifs as limits,^17,24 other studies have concluded that the predominant configurations can be described as either all Zundel-like^26–28 or all Eigen-like.^19,29

Infrared (IR) spectroscopy provides a stringent benchmark to measure molecular configurations^{14,18,23,30–33} owing to its sensitivity to the nuclear potential, but modeling the infrared spectrum of the aqueous proton has proven to be challenging due to the extreme anharmonicities of the nuclear potential.^34–36 This challenge has spurred the development of highly precise quantum mechanical methods such as the multiconfiguration time-dependent Hartree (MCTDH)^34,35 and vibrational self-consistent field/virtual state configuration interaction (VSCF/VCI) methods^36–38 to calculate the IR spectrum of gas-phase protonated water clusters. The highly accurate IR spectra from these techniques have been invaluable for disentangling the behavior of highly mixed anharmonic modes in gas-phase clusters.^{34,35,37,39,40} Due to the expense of these techniques, fully quantum mechanical calculations on high-dimensional potential energy surfaces (PESs) have been limited to small, cold gas-phase clusters H⁺(H₂O)_n, where n < 21.^36,37,40,41

Analysis of the aqueous proton IR absorption spectrum is particularly challenging because it appears as a continuum spanning 1000 cm⁻¹–3000 cm⁻¹.^25,42,43 Molecular dynamics simulations have reproduced the aqueous proton continuum using the Fourier transform of the dipole autocorrelation function,^42,44–46 but this approach has only recently been employed for decomposing the spectrum by Eigen-like or Zundel-like configurations.⁴⁶ To examine the spectral signatures of different aqueous proton configurations, some of us recently analyzed the harmonic IR spectra of thousands of aqueous proton clusters generated from aqueous multi-state empirical valence bond (MS-EVB) simulations.⁴³ While this approach illuminated spectral contributions due to hydration symmetry with reliable statistics, it came at the cost of accuracy from the harmonic approximation.

We recently improved upon those harmonic calculations by employing high-level VSCF/VCI calculations on an ensemble of 800 configurations from MS-EVB simulations.³⁹ To calculate the anharmonic IR spectra of these clusters, we employed the local monomer approximation,⁴⁷ where we chose a “special pair” local monomer H⁺(H₂O)₂ whose spectrum was calculated fully quantum mechanically in an environment of four static solvating water molecules. As a result, we were able to achieve a much more reliable decomposition of the linear IR spectrum and address several controversies in feature assignment. In particular, we were able to reliably distinguish features arising from Eigen-like, Zundel-like, and intermediate configurations in the liquid.

The linear IR spectrum of H⁺(aq) provides a limited amount of information due to the high spectral congestion of the continuum and multiple types of broadening. To address these issues, recent nonlinear IR spectroscopic studies of the aqueous proton have provided a wealth of new data that require further interpretation. The two-dimensional IR (2D IR) spectrum reduces spectral congestion by spreading out the IR response over two frequency axes, distinguishes homogeneous and inhomogeneous spectral broadening mechanisms, and identifies anharmonically coupled pairs of modes in the form of off-diagonal cross peaks.^48–50 The three ultrashort IR pulses used in the experiment provide a 100-fs time resolution and access to higher-lying excited states such as overtones and combination bands.⁴⁸ All of these features make a 2D IR spectrum a highly stringent observable for modeling the nuclear potential and distribution of configurations.

A burst of 2D IR experiments in the last few years has provided key insights into characterizing aqueous proton configurations. For instance, direct evidence of aqueous Zundel-like configurations was established via a cross peak between the flanking water stretching and bending modes at 3200 cm⁻¹ and 1750 cm⁻¹, respectively.⁵¹ The measurement of the 2D IR spectrum along the diagonal at 1200 cm⁻¹ revealed an excited state absorption (ESA) at a higher frequency than the ground state bleach (GSB), providing undeniable evidence of shared protons tightly confined by two flanking waters in Zundel-like configurations.^26,52 Finally, recent 2D IR spectra at multiple key diagonal and off-diagonal regions across the mid-IR displayed cross peaks between all vibrational features, suggesting the coupled vibrations all belong to a single species with a broad distribution of structures in solution, rather than distinct Eigen and Zundel configurations.⁵³ The 2D IR spectrum was most consistent with a distribution of asymmetric Zundel-like vibrational potentials where the excess proton feels the effects of confinement between two waters, but a connection between the potential shape and molecular geometry is missing.

Even with the strengths of 2D IR spectroscopy, the substantial anharmonicity of the nuclear potential has hampered a straightforward analysis of the nonlinear spectroscopic data. For instance, a clear relationship between proton stretch frequencies and geometric coordinates of the aqueous proton complex has not been determined like in water or other aqueous solutions.^30,54 Additionally, the proton stretching vibration of Zundel-like configurations and the umbrella motion of Eigen-like configurations are both resonant at ∼1200 cm⁻¹,^23,39 which introduces ambiguity into the previous interpretation of the 2D IR spectrum. Finally, the polarization relationships between 2D IR features have not yet been connected to relative orientations between vibrational transition dipoles. With an establishment of spectrum–structure relationships, the 2D IR spectrum could be employed as a constraint on possible configurations and inform simulations of the aqueous proton.

To address the above challenges, we have employed anharmonic, high-level VSCF/VCI calculations to analyze the 2D IR spectrum of the aqueous proton. The previous calculations on 800 [H⁺(H₂O)₂](H₂O)₄ clusters include frequency and dipole information for all of the excited states accessed in the 2D IR measurement, so they provide an opportunity to test our previous assignments. In the static limit, the 2D IR spectrum can be equated with a two-dimensional conditional probability frequency distribution function of vibrational frequencies, describing the spectrum of infrared transitions that can be detected, given that a fundamental vibrational transition has been excited. We first present such 2D conditional probability spectral distributions for single configurations of the gas-phase 6-Zundel and 6-Eigen protonated water hexamers, which serve as limiting cases for the distribution of structures seen in the aqueous phase. We then present 2D conditional probability spectra for Zundel-like, Eigen-like, and intermediate aqueous configurations. An analysis of the proton stretch vibrations is presented for interpreting the 2D IR diagonal region around 1200 cm⁻¹, followed by an analysis of the off-diagonal cross peaks due to Zundel-like configurations. These spectral assignments provide a rigorous foundation for structural interpretation of 2D IR spectroscopy of protonated water complexes.

In our previous study on the linear IR spectrum, we carried out VSCF/VCI calculations on 800 aqueous H⁺(H₂O)₆ clusters drawn from the MS-EVB trajectory from Ref. 43. In this work, we used the same VSCF/VCI results to study their excited-state spectra.³⁹ Details pertaining to the MS-EVB simulation⁴³ and VSCF/VCI calculations³⁹ have been described in detail previously, but we briefly present them here. We extracted 800 random frames from a 1-ns trajectory of 256 simple point charge flexible (SPC/Fw) water molecules, one excess proton, and one chloride using the MS-EVB 3.2 model.⁵⁵ From each frame, the excess proton was determined as the H atom on hydronium with the smallest proton sharing parameter δR_OH = |R_O1,H − R_O2,H| on the hydrogen-bonded triplet O₁–H⋯O₂. We previously verified that the δR_OH and R_O1,O2 distributions from the 800 clusters reflected those from the entire trajectory.³⁹ As with our previous study, we acknowledge here that the structural statistics in the MS-EVB model may slightly deviate from experiment, but this simulation is useful for generating a distribution of configurations that one might expect to find in the aqueous phase. For each of the 800 frames, we selected H⁺(H₂O)₆ clusters consisting of the excess proton and its six closest waters for VSCF/VCI calculations. As well, we previously verified that the harmonic spectra of the protonated hexamer clusters closely resembled the spectra of those clusters with up to 18 solvating molecules, indicating that the main spectroscopic features are well represented with hexamer clusters.³⁹ 

VSCF/VCI calculations of the vibrational spectra for each hexamer were calculated on a highly accurate ab initio many-body potential energy surface (PES) and dipole moment surface (DMS), which have been successfully used to calculate the structures, energies, and IR spectra of gas-phase protonated water clusters.^{37,41,56–58} To reduce the computational cost, we employed the local monomer approximation on the 800 aqueous proton configurations, which has previously been successful at modeling the IR spectrum of liquid water and ice.^36,38,47,59 In this scheme, the cluster is partitioned into monomers, and the VSCF/VCI calculation is performed on each monomer while fixing the others, significantly speeding up the computation by reducing the number of coupled normal modes. We used a H⁺(H₂O)₂ local monomer, consisting of a central excess proton and its two closest waters, which was the smallest local monomer that captured both hydronium and Zundel-like geometries. We note that for the purposes of capturing a broad distribution of configurations, we did not relax the local monomer geometry before spectral calculations. Previous vibrational simulations on aqueous proton clusters at the harmonic level demonstrate only a small fraction of modes with imaginary frequencies, which lie below 500 cm⁻¹.⁴³ To generate a basis for describing anharmonic vibrations, the normal modes of the H⁺(H₂O)₂ local monomers were calculated. Using a frequency cutoff of 700 cm⁻¹, ∼12–15 normal modes were retained for each configuration. The calculations were performed with the MULTIMODE software,⁶⁰ using a 4-mode representation of the potential and an excitation space of single, double, triple, and quadruple excitations. From the resulting excited states, we are able to access the frequency positions and transition dipole vectors of fundamental excitations and excited-state transitions into overtones and combination bands that can be compared to the 2D IR experiment.

We begin with calculated conditional frequency probability distributions for the 6-Zundel and 6-Eigen gas-phase hexamers to highlight the impact of anharmonicity on their characteristic spectral features and illustrate the expected 2D IR features associated with these structures. The conditional probability frequency distributions P_Fund(ω₃|ω₁) and P_Exc(ω₃|ω₁) represent the spectrally weighted expectation values of detecting a frequency ω₃ given excitation of a fundamental transition at ω₁^61,62 and are defined here as

These distributions are expressed in terms of angular transition frequency ω_fi associated with the excitation from |i⟩ → | f⟩ with corresponding transition dipole vectors $μ⇀fi$⁠. In Eq. (1a), the summation over a and b refer to singly excited (one-quantum) states, whereas in Eq. (1b) the index a refers to a singly excited state and f can be a singly excited or higher-lying excited state with energy up to 6775 cm⁻¹. P_Exc(ω₃|ω₁) is primarily represented by excitations into two-quantum states such as overtones and combination bands, but the strong anharmonicity of the system allows excitation from one-quantum states to other one-quantum states or higher-lying excited states such as three-quantum states. The intensity of each feature is captured by the product of $μfi2$ terms. Angular brackets denote ensemble averaging over calculations on many static clusters.

The sum of Eqs. (1a) and (1b) is closely related to the static limit of the 2D IR spectrum for these clusters, with no contribution to the line shapes from vibrational dynamics. P_Fund(ω₃|ω₁) represents the positive bleach of fundamental transitions in 2D IR spectra, while P_Exc(ω₃|ω₁) represents the negative excited state absorptions between the first excited state |a⟩ and the final excited state |f⟩. These spectral densities are calculated from the coupled eigenstates from the clusters and are not derived from third-order response function calculations. As such, these surfaces resemble absorptive 2D IR spectra but do not arise from a response function calculation or a sum over Liouville pathways. The set of accessible excited states from a collection of N fundamental transitions includes N(N + 1)/2 overtones and combination bands, which behave more anharmonically and are more strongly mixed than the fundamentals. As a result, the usual harmonic approximation does not hold, i.e., the μ₂₁ and μ₁₀ transition dipoles are not necessarily spatially aligned and the intensity ratio μ₂₁²/μ₁₀² ≠ 2. The large collection of excited states and breakdown of harmonic rules result in many features in P_Exc(ω₃|ω₁) that appear dispersed along ω₃.

The 6-Zundel isomer consists of a symmetric “Zundel” core with four solvating waters. Its absorption spectrum consists of three characteristic peaks in the VSCF/VCI calculated spectrum (Fig. 1): an intense peak at 1102 cm⁻¹ from the asymmetric stretching motion of the central proton (peak 1), a peak at 1765 cm⁻¹ (peak 2) from the concerted bending motion of the two flanking waters in the Zundel core, and an intense and broad doublet peak at 3128 cm⁻¹ and 3207 cm⁻¹ (peak 3) due to symmetric and antisymmetric O–H stretch normal modes of the two waters flanking the central proton, respectively. P_Fund(ω₃|ω₁) for the 6-Zundel cluster [Fig. 2(a)] displays a distinct coupling pattern between these peaks spread across the mid-IR (Table S1). For example, the off-diagonal peaks 4–6 arise from anharmonic coupling between the proton stretch, flanking bend, and flanking stretches, and the analogous downhill cross peaks 7–9 appear across the diagonal. The proton stretch-associated cross peaks are particularly bright, which distinguishes the Zundel spectral pattern.

The 6-Eigen cluster shows a different spectral pattern from the 6-Zundel species. This isomer can be described as an Eigen species H₃O⁺(H₂O)₃ with two waters forming an incomplete second hydration shell around the central hydronium. This arrangement introduces asymmetry in the solvation environment of the hydronium compared to the Eigen cation H₉O₄⁺. Its linear IR spectrum [Fig. 1(d) and Table S2] includes three intense peaks at 1991 cm⁻¹, 2345 cm⁻¹, and 2893 cm⁻¹ (peaks 3′–5′, respectively) from the three local O–H stretches of the central hydronium, in order from strongest to weakest H-bonds. The lowest frequency O–H stretch participates in what is known as the “special pair” of the hydronium and its closest first-shell water. In the 6-E complex, the special pair first-shell water donates a H-bond to a second shell water. The intermediate-frequency O–H stretch involves another branch with first and second shell waters, whereas the high-frequency O–H stretch donates to the first-shell water without a second shell water. The hydronium umbrella and in-plane bending modes appear as weak features at 1209 cm⁻¹ and 1619 cm⁻¹ (peaks 1′ and 2′, respectively). The P_Fund(ω₃|ω₁) distribution [Fig. 2(b)] reveals a high cross peak density between 2000 cm⁻¹ and 3300 cm⁻¹ along both axes [inside the dashed box in Fig. 2(b)].

The P_Exc(ω₃|ω₁) distributions demonstrate that excited states are accessible from the fundamental transitions in the 6-Zundel and 6-Eigen clusters. For the 6-Zundel cluster, the main peaks are associated with overtones or combination bands that arise from the couplings between the fundamental modes in the 6-Zundel cluster. As a result, P_Exc(ω₃|ω₁) displays a similar pattern to P_Fund(ω₃|ω₁) but with frequency shifts due to the anharmonicity of the nuclear potential. While P_Fund(ω₃|ω₁) is symmetric across the diagonal, P_Exc(ω₃|ω₁) necessarily loses that symmetry since excitations into two-quantum states are dependent on the initially accessed one-quantum state. Most anharmonic ESAs are lower in frequency than the GSB,^48,63 but the ESAs of all of the 6-Zundel vibrations appear above the diagonal. This is readily observed in the combined P_Fund(ω₃|ω₁) + P_Exc(ω₃|ω₁) surface in Fig. 2(e), which is the closest available match to the 2D IR spectrum and includes the interference of positive and negative features. Combination bands between the proton stretch, flanking bends, and flanking stretches are accessed off the diagonals as peaks 4^*–9^*. In the 6-Eigen cluster [Fig. 2(d)], there is a high density of excited states between 2000 cm⁻¹ and 3000 cm⁻¹ and forms a similar box pattern like the fundamentals in Fig. 2(b). However, assignments of these features are challenging because there is strong anharmonic mixing between overtones and combination bands. For example, the ESA at (ω₁, ω₃) = (2000 cm⁻¹, 2080 cm⁻¹) arises from the excitation from the special pair O–H stretch to a combination band involving the special pair O–H stretch overtone and another combination band between the hydronium umbrella mode and hydronium libration.

There are some features about the P_Exc(ω₃|ω₁) distributions that are less intuitive than for the fundamental distributions due to the large anharmonicities in these clusters. For both clusters, the intensities of the features in P_Exc(ω₃|ω₁) differ from their counterparts across the diagonal, which is one indication of electrical anharmonicity in these clusters and the large manifold of the accessible excited states. Here, excitation of the proton stretch amplifies the intensities of other excited-state transitions, compared to the reverse excitation order. Many ESA peaks also appear elongated over hundreds of wavenumbers along ω₃ because of excitations into various combination bands and overtones. The proton stretch and flanking bend modes in the 6-Zundel cluster display substantial anharmonic character: the μ₂₁ transition dipoles are 9° and 7° misaligned, respectively, from the corresponding μ₁₀ vectors, and the intensity ratios μ₂₁²/μ₁₀² are 0.286 and 0.310, respectively, instead of the expected value of two for harmonic transition dipoles. Finally, harmonically forbidden direct excitations from the proton stretch and flanking bend into the flanking stretch fundamental are present in Fig. 2(c) as peaks A and B, respectively. P_Exc(ω₃|ω₁) surfaces with only two-quantum transitions are provided in Fig. S1.

To study the hydrated proton clusters in more detail, we next analyze the vibrational frequencies assigned to the asymmetric proton stretch in 6-Z and special pair O–H stretch in 6-E since the proton stretch has been shown to shift dramatically as a function of hydration in gas-phase clusters.^32,64 As in prior work, we describe the vibrations in each cluster in a harmonic basis, assigning the proton stretch normal mode to the mode that showed the most displacement of the local monomer central proton. The anharmonic proton stretch eigenstate was selected as the vibration with the largest squared VCI coefficient in the proton stretch normal mode coordinate.³⁹ The anharmonicity pattern of the fundamental (ω₁₀) and excited state (ω₁₂) transitions can be used to categorize the nuclear potential. In a harmonic system, the frequency ratio γ = ω₂₁/ω₁₀ = 1.0 arises from the parabolic nuclear potential. We therefore classify anharmonic potentials as superharmonic if the transition frequency ratio γ = ω₂₁/ω₁₀ > 1 or subharmonic if γ < 1. An example of a superharmonic potential shape is the quartic potential in one coordinate V(r) = χr⁴, where γ = 1.33. γ is a similar measure to the anharmonicity parameter Δ = ω₂₁ − ω₁₀ but is inherently referenced to the magnitude of ω₁₀. While assigning a potential to the single number γ ignores many details, it is practical because it is directly connected to experimental observables and can be calculated from an arbitrary potential shape.

The ω₂₁ > ω₁₀ anharmonicity pattern [Fig. 3(a)] of the 6-Zundel proton stretch provides indication of its superharmonic potential shape. For illustrative purposes, we show the slice of the 6-Zundel nuclear potential along the proton stretch normal mode coordinate in Fig. 3(b). The potential shape consists of a relatively flat bottom with steep, confining walls on both sides compared to the harmonic oscillator, which raises the energy of the second excited state. For the 6-Zundel proton stretch, γ = 1.324. By contrast, the ESA frequency of the 6-Eigen special pair O–H stretch at 1991 cm⁻¹ is slightly less than its respective fundamental transitions [Figs. 3(c) and 3(d)] such that γ < 1. Interestingly, there is no “pure” O–H stretch overtone in the 6-Eigen hexamer because the overtones are highly mixed with other excited states. There are also multiple combination bands that mix with the overtone, as seen above the 6-Eigen proton stretch fundamental [Fig. 3(c)]. Importantly, the proton stretch potentials shown in Fig. 3 are one-dimensional slices on a multidimensional potential energy surface that couples the proton stretch to other degrees of freedom such as bending and O⋯O extension. In Fig. S2, we include two-dimensional slices of the nuclear potential as a function of proton stretch, flanking bending, and O⋯O extension for the 6-Z and 6-E clusters.

We previously analyzed the proton stretch 2D frequency spacing and line shape in terms of an asymmetric Zundel-like quartic potential shape in one dimension in terms of the proton stretch normal mode coordinate r,⁵³ 

We found that the quartic confinement coefficient χ and the linear asymmetry β trend strongly with proton stretch frequencies, whereas the quadratic barrier term κ is small compared to the proton stretch zero-point energy. The magnitude of χ reflects the confinement of the shared proton between two flanking waters in the 6-Zundel hexamer. As a result, the potential shape pushes up the second excited-state energy. Increasing β drives up ω₁₀ (as seen between 6-Zundel and 6-Eigen), while ω₂₁ is not as strongly affected, resulting in a lower γ value in the 6-Eigen complex.

Orientational information between two transition dipoles $μ⇀1$ and $μ⇀2$ can be obtained from polarization-resolved experiments by measuring the 2D anisotropy $r(ω1,ω3)=0.4P2μ^1(ω1)⋅μ^2(ω3)$⁠, where P₂ is the second Legendre polynomial, circumflexes over the dipoles indicate unit vectors, and the angular brackets denote ensemble averaging in isotropic media. For two parallel transition dipoles, r = 0.4, whereas two perpendicular dipoles produce an anisotropy of r = −0.2. Here, we calculated the polarization-resolved parallel and perpendicular surfaces similarly to the P_Fund(ω₃|ω₁) and P_Exc(ω₃|ω₁) surfaces in Eqs. (1a) and (1b), including orientational correlation prefactors Y_∥,⊥(θ),

where $Y∥(θ)=1151+2⁡cos2⁡θ,$$Y⊥(θ)=1152−cos2⁡θ,$ and $θ=μ^1(ω1)×μ^2(ω3)$⁠.⁶⁵ With polarization-resolved surfaces P_∥ and P_⊥, the 2D anisotropy can be calculated from the parallel and perpendicular probability surfaces as

The 2D anisotropy can be calculated for P_Fund(ω₃|ω₁), P_Exc(ω₃|ω₁), or the combined surface P_Fund(ω₃|ω₁) + P_Exc(ω₃|ω₁). For the combined surface, it is important to calculate the combined P_Fund,∥(ω₃|ω₁) + P_Exc,∥(ω₃|ω₁) and P_Fund,⊥(ω₃|ω₁) + P_Exc,⊥(ω₃|ω₁) surfaces separately before calculating anisotropy in order to properly account for interference between positive and negative features.

The 2D anisotropy surfaces for the 6-Zundel and 6-Eigen clusters display distinct anisotropy patterns that further distinguish these two clusters. In the 6-Zundel cluster [Fig. 4(a)], the anisotropies of diagonal peaks 1 and 2 are both 0.4 since the modes along the diagonal are parallel with themselves. However, the anisotropy of peak 3 is less than 0.4 likely due to interference of anisotropies between symmetric and antisymmetric flanking stretches. The 2D anisotropies of the cross peaks encode the relative orientations between two vibrational modes. For instance, the anisotropy of peak 4 is near 0.4, indicating that strongly coupled proton stretch and flanking bend modes are roughly parallel to each other, as expected from symmetry considerations. Peak 5 shows regions of both positive and negative anisotropy, indicating that the proton stretch is parallel to the symmetric stretch (3128 cm⁻¹) and perpendicular to the asymmetric flanking O–H stretch (3207 cm⁻¹). The 2D anisotropy plot for the 6-Eigen complex [Fig. 4(b)] demonstrates highly varying orientational information between the hydronium O–H stretches. Along the diagonal, each stretch displays an anisotropy of 0.4, but the anisotropies of their cross peaks are close to zero or negative, reflecting the large angles between them. On the other hand, the stretch-umbrella cross peaks all show anisotropy of −0.2, which is consistent with the out-of-plane transition dipole of the umbrella bend. The 2D anisotropy patterns of the 6-Zundel and 6-Eigen P_Exc(ω₃|ω₁) and combined 2D anisotropy surfaces [Figs. 4(c)–4(f)] are also similar to those for the Zundel and Eigen fundamental surfaces.

To analyze the conditional probability frequency distributions of aqueous proton configurations, we calculated the anharmonic IR spectral densities for the H⁺(H₂O)₂ local monomers within 800 H⁺(H₂O)₆ clusters drawn from MS-EVB configurations and analyzed them relating geometric and spectral characteristics. We previously found that ⟨R_OH⟩, the quantum expectation value of the distance between the shared excess proton and its nearest oxygen, strongly correlated with the stretch frequency of the excess proton,³⁹ whereas grouping by ⟨R_OO⟩ did not differentiate the spectra as well (Fig. S3). We also found that grouping clusters based on ⟨R_OH⟩ produced distinct spectra that resembled the 6-Zundel and 6-Eigen gas-phase clusters in the extreme limits. In this vein, we analyze structurally averaged spectra by grouping clusters into three groups: Zundel-like (⟨R_OH⟩ > 1.14) Å, Eigen-like (⟨R_OH⟩ < 1.07 Å), and intermediate (1.07 Å < ⟨R_OH⟩ < 1.14 Å). These cutoffs represent one standard deviation away from the median ⟨R_OH⟩ value (1.11 Å), and different choices of cutoffs do not qualitatively change the analysis. In our previous work, we found a roughly linear correlation between ⟨R_OH⟩ and the proton sharing parameter ⟨δR_OH⟩.³⁹ For reference, ⟨R_OH⟩ > 1.14 Å roughly corresponds to ⟨δR_OH⟩ < 0.15 Å, and ⟨R_OH⟩ < 1.07 Å roughly corresponds to ⟨δR_OH⟩ > 0.30, consistent with previous definitions of Zundel-like and Eigen-like geometries.^43,66 For each of these configurations, distributions of ⟨R_OO⟩ distances are broad and significantly overlap (Fig. S4), but the median values only slightly shift from 2.38 Å, 2.41 Å, to 2.43 Å for Zundel-like, intermediate, and Eigen-like configurations, respectively.

For Zundel-like configurations, the linear spectral density (Fig. 5) resembles that of the 6-Zundel gas-phase cluster [Fig. 1(a)] with three features peaked at 1300 cm⁻¹, 1800 cm⁻¹, and 3100 cm⁻¹ and corresponds to the proton stretch, flanking water bends, and flanking O–H stretches in aqueous configurations, respectively.^39,43,51,67 On the other end of the distribution, the aqueous Eigen-like spectrum has features that resemble the gas-phase 6-Eigen spectrum: a weak umbrella mode at 1300 cm⁻¹ and a broad band of hydronium and first-shell water O–H stretches spanning 2000 cm⁻¹–3500 cm⁻¹. The spectrum of intermediate configurations displays three bands like the Zundel-like configurations but is substantially broadened and contributes to the continuum of absorption between 2000 cm⁻¹ and 3000 cm⁻¹ like the Eigen-like geometries. The relative populations of Zundel-like and Eigen-like configurations sampled may not be exactly representative of the bulk,³⁹ so spectral densities and two-dimensional conditional probability surfaces are all presented, normalized by the number of clusters in each geometric grouping. Between the three groups, the integrated intensity of the normalized Zundel-like spectral density is highest, whereas the intermediate spectral density is ∼96% as intense and Eigen-like configurations are ∼85% as intense. Because the Eigen-like features are spread over one broad frequency range instead of appearing in three peaks, the maximum intensity is somewhat weaker than the maximum in the Zundel-like spectral density.

In Fig. 6(a), we present the conditional probability frequency distribution of fundamental transitions in primarily Zundel-like clusters. We observe a diagonally symmetric distribution with features that retain the character of the peaks in the 6-Zundel P_Fund(ω₃|ω₁) surface [Fig. 2(a)]. The diagonal region between 1100 cm⁻¹ and 1500 cm⁻¹ arises from the fundamental transitions of the proton stretching modes and is the most intense feature of the conditional frequency distribution. There is also a clear cross peak pattern between the proton stretch, flanking bend, and flanking stretch (Table S4). The intermediate P_Fund(ω₃|ω₁) [Fig. 6(b)] displays many of the trends of the Zundel-like distribution, but the line shapes are generally broader. On the other extreme, the Eigen-like complexes show a P_Fund(ω₃|ω₁) pattern [Fig. 6(e)] that forms a box between 2000 cm⁻¹ and 3000 cm⁻¹ like the gas-phase 6-Eigen complex (Table S5).

The ESA conditional frequency distributions P_Exc(ω₃|ω₁) are presented in Figs. 6(d)–6(f). For Zundel-like configurations [Fig. 6(d)], the stretch and bend ESAs (peaks 1^* and 2^*, respectively) appear above the diagonal like in the 6-Zundel cluster, whereas the flanking O–H stretch ESAs (peak 3^*) appear along the diagonal and red-shifted compared to the gas phase. The proton stretch ESA and proton stretch/flanking bend combination band (peak 4^*) both display diagonal elongation, which indicates inhomogeneous broadening in their line shapes. The intermediate ESA P_Exc(ω₃|ω₁) retains the same features [Fig. 6(e)], but the overtone near 2000 cm⁻¹ is more intense due to the presence of blue-shifted proton stretch vibrations. Two-quantum transitions are also present in these surfaces but are less distinct than in the gas phase [Fig. 2(c)] and are shown separately in Fig. S5. The Eigen-like excited-state spectral density [Fig. 6(f)] displays the same high spectral density at (ω₁, ω₃) = (2000 cm⁻¹–3000 cm⁻¹, 2000 cm⁻¹–3000 cm⁻¹) like in the P_Fund(ω₃|ω₁) for Eigen-like configurations. As a result, there is a significant destructive overlap in the combined surface for Eigen-like clusters, resulting in low probability amplitude across both dimensions in the mid-IR [Fig. 6(i)].

The 2D anisotropies (Fig. 7) in these groupings mirror those of the gas-phase clusters, which reinforces their separation by Zundel-like and Eigen-like character. For Zundel-like and intermediate configurations [Figs. 7(a) and 7(b)], peaks 4 and 7 show high positive anisotropy (r ≈ 0.25), while the cross peaks to the flanking stretches (peaks 5, 6, 8, and 9) display a primarily negative anisotropy. In the Eigen-like sub-ensemble, the 2D anisotropy pattern varies significantly between 2000 cm⁻¹ and 3000 cm⁻¹, but the cross peaks between the O–H stretches and the umbrella modes are consistently negative. For all three types of configurations, the ESA 2D anisotropies [Figs. 7(d)–7(f)] all show similar trends as in the fundamental 2D anisotropy surfaces. For the combined P_Fund(ω₃|ω₁) + P_Exc(ω₃|ω₁) surfaces, the qualitative trends of positive and negative anisotropy are similar to the individual surfaces, but the effect of interference between positive and negative features results in actual anisotropy values beyond the usual limits of 0.4 and −0.2. To examine the connection between structures and spectra in more detail, we proceed to analyze the piece of the 2D IR spectrum pumped at 1200 cm⁻¹.

Analyzing the 2D diagonal region between 1000 cm⁻¹ and 1600 cm⁻¹ provides an opportunity to learn about the ensemble of Zundel-like and intermediate aqueous proton configurations. To demonstrate this, we present the experimental 2D IR spectrum of this region in Fig. 8(a) from Ref. 53. The ground state bleach lies along the diagonal, while the |1⟩ → |2⟩ ESA appears above the diagonal, centered at ω₃ = 1360 cm⁻¹. This doublet pattern demonstrates a peculiar anharmonicity but is consistent with its assignment as the proton stretch in Zundel-like species.^26,53 In Figs. 8(b)–8(d), we find that the calculated P_Exc(ω₃|ω₁) distributions in this region for Zundel-like, intermediate, and Eigen-like clusters [Figs. 8(b)–8(d)] share an ESA above the diagonal and also display diagonal elongation. For Zundel-like and intermediate clusters, the ESA is primarily due to the proton stretch ω₂₁ transition, whereas the ESA in Eigen-like clusters [Fig. 8(d)] arises from the ω₂₁ transition of hydronium umbrella mode, which has a much smaller transition dipole amplitude. When normalized by the number of clusters in each type of configuration, the Zundel-like and intermediate surfaces are on average significantly more intense than the Eigen-like clusters in this region.

The diagonally elongated probability densities of Figs. 8(b) and 8(c) indicate that there are frequency correlations within a large distribution of ω₁₀ and ω₂₁ transitions. To connect proton stretch frequencies to atomic displacements, we employed linear least-squares regression on ω₁₀ and ω₂₁ with collective geometric coordinates.⁶⁸ We began by optimizing the R² correlation coefficient between the vibrational frequency and a mixed coordinate q modeled by a linear combination of quantum expectation values of local monomer atomic distances, H–O–H angles, and H–O–O–H dihedrals, allowing their weighting coefficients to float freely. We found that the two most significant parameters were the ⟨R_OH⟩ and ⟨R_OO⟩ distances and therefore reoptimized the correlations on a reduced collective coordinate,

where the absolute values of c₁ and c₂ were constrained to sum to one. We carried out this analysis separately for Zundel-like, Eigen-like, and intermediate cluster geometries.

Table I shows the retrieved c₁ and c₂ coefficients and the R² value for corresponding scatter plots in Fig. 9. We find that in Zundel-like configurations ⟨R_OO⟩ is a more influential predictor of ω₁₀ than ⟨R_OH⟩, with |c₁| = 0.40 and |c₂| = 0.60, where both shorter ⟨R_OH⟩ and ⟨R_OO⟩ distances drive up the proton stretch frequency. At the other extreme, the frequency trend in Eigen-like clusters is dominated by ⟨R_OH⟩ (c₁ = 0.96). As seen in Table I and Fig. 9(d), c₁ and c₂ values are similar for ω₁₀ and ω₂₁ trends in each grouping, which points to frequency correlation between ω₁₀ and ω₂₁ and diagonal elongation in the surfaces in Fig. 8. Interestingly, the frequency trends for the intermediate complexes are dominated by ⟨R_OH⟩, even though the spectral density of intermediate complexes more closely resembles that of the Zundel-like configurations (Fig. 5). This suggests that the intermediate complexes represent a turning point between two regimes dominated by ⟨R_OH⟩ and ⟨R_OO⟩. The influence of ⟨R_OO⟩ on the proton stretch frequency in Zundel-like and intermediate configurations indicates that confinement from the flanking two waters significantly influences the shape of the proton stretch potential. However, because of the substantial overlap between proton stretch frequency distributions for Zundel-like and intermediate configurations, the broad distribution of ω₁₀ and ω₂₁ frequencies between ω₁ = 1000 cm⁻¹ and 1600 cm⁻¹ reflects the spread of both ⟨R_OO⟩ and ⟨R_OH⟩ in Zundel-like and intermediate configurations.

While the ω₁ = 1000 cm⁻¹–1600 cm⁻¹ region primarily reports on the heterogeneity of Zundel-like and intermediate distributions, trends across the entire mid-IR can be used as an experimental handle to characterize the heterogeneity of Eigen-like configurations. Figure 9(c) demonstrates that proton stretch ω₁₀ spans almost 1000 cm⁻¹ due to the ⟨R_OH⟩ range of ∼0.06 Å in Eigen-like configurations. This is a relatively small distribution of values, particularly compared to the typical values of ⟨δR_OH⟩ (0.0 Å–0.5 Å), underscoring the extreme sensitivity of the O–H stretch on the bond length. Despite the extreme red-shifting of the special pair O–H stretch, the off-diagonal intensities in the P_Fund(ω₃|ω₁) and P_Exc(ω₃|ω₁) surfaces of Eigen-like clusters [Figs. 6(c) and 6(f), respectively] demonstrate that all three of hydronium’s O–H stretches are anharmonically coupled to each other. This implies that while the special pair H-bond is the primary determinant of the special pair O–H frequency, the rest of the solvation environment will have some influence on the special pair frequency as well.

As discussed in Sec. IV B, the proton stretch frequency is a highly sensitive reporter of the configuration of the aqueous proton complex. In the static limit presented here, ω₁₀ spans ∼2000 cm⁻¹, spectrally separating configurations with different ⟨R_OH⟩ and ⟨R_OO⟩ values. We can compare the ω₁₀ and ω₂₁ values for each cluster to uncover trends in the proton stretch potential energy surface. To start, we present a scatter plot comparing proton stretch ω₁₀ frequencies with their corresponding ω₂₁ transitions (Fig. 10). As shown in Sec. IV B, the ω₁₀ transitions of Zundel-like and intermediate clusters overlap significantly and are distributed across 900 cm⁻¹–1600 cm⁻¹, whereas Eigen-like special pair O–H stretches are almost all above ω₁₀ > 1600 cm⁻¹. The ω₂₁ transitions for the Zundel-like and intermediate clusters are almost all at higher frequency than their corresponding ω₁₀ transitions, placing them above the diagonal in this plot. This is consistent with the surfaces shown in Figs. 8(b) and 8(c) for these two types of clusters. By averaging the ω₂₁ values over 50 cm⁻¹ intervals in ω₁₀, we can calculate a mean frequency ratio ⟨γ⟩ = ⟨ω₂₁⟩/⟨ω₁₀⟩. In Fig. 10, ⟨γ⟩ ≈ 1.3 when ω₁₀ < 1000 cm⁻¹ and gradually reduces to ⟨γ⟩ ≈ 1.15 for ω₁₀ ≈ 1200 cm⁻¹, in a range consistent with the experimental surface in Fig. 8(a).⁵³ For Eigen-like clusters, the ω₂₁ transitions are closer to their corresponding ω₁₀ transitions approaching ⟨γ⟩ ≈ 1.

The above trends in frequency spacing derive from the shapes of the nuclear potential. Figure 11 shows representative one-dimensional potential energy slices along the proton stretch normal mode coordinate defined in Sec. III B. The potential energy curve of the Zundel-like configuration [Fig. 11(a)] resembles the 6-Zundel potential energy curve in Fig. 3(b), featuring a relatively flat bottom with steep walls and a mostly symmetric shape. These are the characteristics of a superharmonic potential shape, leading to the ω₂₁ > ω₁₀ (γ > 1) relationship seen in Fig. 10. The potential of the intermediate complex [Fig. 11(b)] shows similar steep walls, but the bottom of the potential appears more asymmetric. Following the discussion of the relative effects of confinement and asymmetry in Sec. III B, the increase in asymmetry pushes up the ω₁₀ transition, reducing γ to about 1.15. At the high-frequency extreme, the Eigen-like potential [Fig. 11(c)] displays a highly asymmetric shape much like the 6-Eigen potential [Fig. 3(d)] resulting in γ < 1.

While these potentials are representative for Zundel-like, Eigen-like, and intermediate configurations, there is some variation within each category (Fig. S6). For each geometric group, one can find examples of flat-bottomed potentials, potentials with substantial asymmetry, and there are even a couple of cases of low-barrier double-welled potentials in each group. There are only a handful of such potentials in our ensemble, whereas the vast majority of complexes show single-well potential shapes along the proton stretch normal mode coordinate.

Each of the presented potentials are only slices in one dimension, but the full potential is multidimensional and involves coupling to other degrees of freedom (Fig. S7). VSCF/VCI calculations reveal that the quantitative calculation of proton stretch frequency requires at least four dimensions, whereas low-dimensional surfaces fail to reproduce frequency positions.^36,38,39 Anharmonic coupling with other modes influences the observed ω₁₀, ω₂₁, and γ in non-trivial ways, which could explain the difficulty in achieving strong correlations between the proton stretch and geometric descriptors such as ⟨R_OH⟩ and ⟨R_OO⟩. Strong anharmonic mixing with other modes made identifying the proton stretch ambiguous in some cases. These issues could be alleviated by employing a similar analysis on a simpler strongly H-bonded aqueous system with clearly delineated vibrational modes. As a result, linear regression analyses correlating collective geometric coordinates with proton stretch potential characteristics like confinement or asymmetry did not yield quantitatively satisfactory correlations, even though the trends based on Eq. (2) were physically intuitive. Finally, the wide variability in one-dimensional potential energy cuts (Fig. S6) is amplified in two-dimensional surfaces, which complicates any systematic analysis of potential energy shape in multiple coordinates. This implies that the parameter γ is influenced by more than the potential shape in one dimension. Even though γ is a conceptually useful parameter, it cannot be connected quantitatively to the potential shape.

There has been illuminating work demonstrating that the proton stretch is a function of many strongly coupled coordinates such as O⋯O stretching, flanking in-plane bending, and proton-involved low-frequency modes in Zundel-like and Eigen-like complexes.^34,35 The VSCF/VCI calculations are able to capture both of these aspects, but the accuracy is only achieved by applying quantum mechanical calculations on high-dimensional potential energy surfaces, each of which can be highly nonintuitive.³⁸ Zundel-like complexes exhibit strong mode mixing in particular, perhaps an effect of the smaller ⟨R_OO⟩ distances which would increase interactions between molecules in the complex. As configurations become more Eigen-like, the proton stretch frequency trends more linearly with the ⟨R_OH⟩ length. This suggests that coupling to other types of modes is less important to predict frequency, but as discussed in Sec. IV B, there is still strong coupling between the three hydronium O–H oscillators.³⁶ 

The conditional frequency distributions presented in Fig. 6 display distinct cross peak patterns that distinguish spectral response from Eigen-like, Zundel-like, and intermediate geometries. For aqueous Zundel-like solutions, the off-diagonal features are relatively intense and tightly distributed in the four corners of the combined P_Fund(ω₃|ω₁) + P_Exc(ω₃|ω₁) surface [Fig. 6(g)]. Here, we examine two of these regions in more detail to interpret the response from experimental 2D IR spectra.

In Fig. 12(a), the isotropic 2D IR spectrum of 2M HCl–H₂O at (ω₁, ω₃) = (1200 cm⁻¹, 1750 cm⁻¹) consists of a broad bleach spanning more than 250 cm⁻¹ in the excitation dimension and ∼150 cm⁻¹ along the detection dimension. Additionally, there are narrow positive and negative bands below ω₃ = 1700 cm⁻¹ due to the anharmonic coupling between high-frequency librational modes and bending modes in bulk water.⁶⁹ The 2D IR spectrum of H₂O in this region is provided in Fig. S8. The bleach at ω₃ = 1750 cm⁻¹ is joined by a negative feature centered at ω₃ = 1900 cm⁻¹. The combined conditional probability surfaces from VSCF/VCI calculations [Fig. 12(b)] indicate that the most intense response originates from Zundel-like geometries with some contribution from intermediate configurations. This comparison confirms our previous assignment⁵³ that the cross peak doublet arises from coupling between the shared proton stretch and concerted flanking bend in Zundel-like configurations.

Closer inspection of these surfaces reveals orientational and frequency correlations between the proton stretch and flanking bend modes. The parallel polarization preference in the experimental spectrum⁵³ and calculated 2D anisotropy surfaces (Fig. S9) reveal that these two modes have roughly parallel transition dipoles, implying that the concerted flanking bend dipole lies along the shared O⋯O axis and has asymmetric normal mode character. Additionally, the calculated conditional frequency probability surface [Fig. 12(b)] displays diagonal elongation of bleach and induced absorptions, indicative of correlated frequencies with the excited proton stretching mode. Inspection of the individual P_Fund(ω₃|ω₁) and P_Exc(ω₃|ω₁) surfaces in Figs. 6(a) and 6(d) reveals that while there is some diagonal elongation of the bleach feature, the majority of the diagonal elongation is due to the frequency correlation in the two ESA features flanking the bleach. This frequency correlation is not as strongly present in the experimental spectrum, however, which suggests fast frequency fluctuations within 150 fs.

Another region of interest is in the cross peak in the (ω₁, ω₃) = (3200 cm⁻¹, 1750 cm⁻¹) region [Fig. 12(c)], which was previously measured as the first experimental evidence for the Zundel species in solution.⁵¹ This cross peak is reproduced in the conditional probability distribution for Zundel-like configurations [Fig. 12(d)], reinforcing this assignment as a cross peak between the Zundel flanking stretches and flanking bends. The negative features flanking the bleach in the conditional probability distribution are due to excitations from the flanking stretch into the flanking stretch/flanking bend and flanking stretch/shared proton stretch combination bands. These are also expected in the experimental 2D IR spectrum, but this region also contains an intense negative background due to the strongly anharmonic water O–H stretch ESA that spans ω₃ = 1500 cm⁻¹–3000 cm⁻¹.⁷⁰ 

In contrast to the clear cross peaks of Zundel-like configurations, the 2D IR features from Eigen-like geometries are less noticeable. Inspecting the P_Fund(ω₃|ω₁) and P_Exc(ω₃|ω₁) surfaces separately [Figs. 6(c) and 6(f), respectively] suggests that there are broad features due to coupling with the Eigen-like stretch continuum, but the combined surface [Fig. 6(i)] reveals that these features destructively overlap. For example, if one examined the P_Fund(ω₃|ω₁) surface alone, it would appear that the cross peak between the Eigen umbrella and O–H stretch modes should be apparent at (ω₁, ω₃) = (1200 cm⁻¹, 2000cm⁻¹–2800 cm⁻¹), but this cross peak disappears when overlapped with the P_Exc(ω₃|ω₁) surface. This is inconsistent with the presence of a cross peak bleach in the experimental 2D IR spectrum in this region.⁵³ One possibility is that the bleach and ESA overlap and destructively interfere in experiment. This would imply that the bleach is overall more intense than the induced absorption, which could be justified by the reduced ESA transition dipole or larger breadth of allowable ESA transitions due to the anharmonicity of the system.

The difference between the calculated and experimental surfaces in this region means that there are a few possible assignments for this cross peak. Previously, some of us measured a cross peak between the 1200 cm⁻¹ and 2500 cm⁻¹ features with perpendicular polarization preference, which was used to assign the 2000 cm⁻¹–3000 cm⁻¹ region as flanking O–H stretches of Zundel-like species.⁵³ However, the calculated linear spectral density and two-dimensional conditional frequency probability distribution point to this region belonging to Eigen-like geometries. Assuming static configurations in solution, this would suggest that there should not be a cross peak between Eigen-like and Zundel-like geometries.

The VSCF/VCI calculations on clusters pulled from MS-EVB snapshots capture the broad distributions of structural parameters like ⟨R_OH⟩ and ⟨R_OO⟩ but do not account for solution fluctuation dynamics. These distributions arise from thermal fluctuations at room temperature,⁵² the low or non-existent barrier for the proton to rattle between water molecules,^24,25 and the quantum delocalization of the excess proton with finite zero-point energy.²⁰ Even in classical MS-EVB simulations, Eigen-like configurations are stabilized compared to Zundel-like configurations by only ∼0.8 kcal/mol,⁵⁵ which means that the excess proton can experience a wide distribution of local environments at room temperature. This breadth in configurations underscores the importance of considering the ensemble of aqueous proton structures as a distribution, rather than characterizing it only as the Eigen and Zundel species that represent the limits of this distribution.

Within this broad distribution, simulations suggest that there are fast geometric fluctuations on the sub-100 fs timescale. Multiple simulations have found sub-100 fs proton position fluctuations,^26,52 including phenomena known as proton rattling⁷¹ and the special pair dance.²⁹ These both impact the instantaneous values of ⟨R_OH⟩ and ⟨δR_OH⟩, which have been used to characterize complexes as Eigen-like or Zundel-like. Recent ab initio and quantum mechanics/molecular mechanics (QM/MM) simulations of the excess proton found that δR_OH quickly decorrelates on a sub-100 fs timescale,^45,52 which suggests that ⟨R_OH⟩ also decorrelates rapidly. This implies that while ⟨R_OH⟩ is spectrally sensitive to instantaneous configuration, this vibration likely will not be sensitive to dynamics associated with slower processes like irreversible proton transfer, which occurs on a 1 ps–2 ps timescale.^72,73 Additionally, the fast and large-amplitude fluctuation of ⟨R_OH⟩ implies that homogeneous broadening is likely a significant component of the linear and nonlinear spectra that cannot be captured with the current VSCF/VCI calculations. Indeed, this is the reason that diagonal GSB features in our calculated spectra, which are by definition perfectly correlated, are not useful in comparing to experiments with massively broadened spectral features. This plays a particularly large role in interpreting Eigen-like clusters in the 2000 cm⁻¹–3000 cm⁻¹ window.

These rapid nuclear fluctuations of the excess proton in an essentially barrier-free potential would imply that the standard adiabatic assumptions for interpreting the spectroscopy are not valid. For instance, the sub-100 fs structural interconversions would occur on a timescale comparable to the vibrational period of the proton stretch mode (∼25 fs). This would imply that the electronic and nuclear degrees of freedom that determine anharmonic coupling are also evolving on the timescale of the vibrational period. This could play an important role in shaping cross peaks between Eigen-like, indermediate, and Zundel-like spectral features at early waiting times. These vibrational non-adiabatic dynamics are known to efficiently dissipate energy in the aqueous medium through collective low-frequency modes,^74,75 and these effects are sure to influence the environment around the excess proton. Sub-Å positional fluctuations of the solvent atoms, which usually are accounted for in models for vibrational line shape broadening, can induce large configurational changes on the characteristic length scale of the excess proton. These factors imply that static configurations can be misleading as true representations of the aqueous system, despite the access to nuclear anharmonicities.

Because of the essentially zero-barrier proton fluctuations, there may be more appropriate geometric coordinates to consider when addressing the interconversion between Zundel-like and Eigen-like geometries. For instance, R_OO has been used to discriminate Eigen-like and Zundel-like configurations in simulations,²⁴ and this coordinate is expected to fluctuate on a 100-fs timescale, as suggested by time-domain THz results.⁷⁶ Even then, if Eigen-like and Zundel-like configurations interconvert on this timescale, there may be a configurational correlation on a ps timescale due to the constraints imposed on the special pair complex by the H-bonding pattern of the solvation environment. Some of us have measured a 2.5-ps anisotropy decay timescale for the Zundel-like bending modes due to Grotthuss proton transfer,^73,77 which is intimately tied to the collective H-bond reorganization process of liquid water.^72,73 These data suggest that fs-timescale thermal fluctuations not totally scramble configurations of the proton hydration complex, even if these fluctuations cause large-amplitude excursions of the excess proton itself.

One limitation of analyzing an ensemble of static configurations is the inability to characterize the excess proton ensemble in a dynamical sense. From one perspective, some studies advocate that the aqueous proton complex is fundamentally Zundel-like, and deviations from the symmetric potential shape come from fast, large amplitude electric field fluctuations.^12,25,26 On the other hand, many studies continue to describe the aqueous proton as a fundamentally hydronium-like or Eigen-like species with large fluctuations in geometry due to solvation dynamics.^{29,71,76,78,79} Each of these conceptions of the aqueous proton can produce the broad distribution of ⟨R_OH⟩ and ⟨R_OO⟩ values that we observed in our study, and perhaps the most appropriate description so far comes from early simulation studies describing it as a “fluxional complex” or dynamic, strongly bonded and persistent “special pair.”²⁴ 

We have presented an analysis of two-dimensional conditional frequency probability distributions generated from anharmonic vibrational calculations on gas-phase H⁺(H₂O)₆ clusters and 800 aqueous proton configurations, which can be used to interpret the 2D IR spectrum of the aqueous proton in the static limit. The gas-phase clusters provide single-structure model systems to interpret how their anharmonic potential energy surface is reflected in two-dimensional conditional frequency probability distributions for transitions between 0 cm⁻¹ and 6800 cm⁻¹. By grouping clusters based on Zundel-like, Eigen-like, and intermediate geometric sub-ensembles according to ⟨R_OH⟩ values, we uncover trends in vibrational frequencies, relative orientations, and potential shape. The vibrational modes of the aqueous proton complex are sensitive reporters of instantaneous geometric configuration, particularly the expectation value of the excess proton-oxygen distance ⟨R_OH⟩ and ⟨R_OO⟩ for Zundel-like configurations.

Analysis of the P_Fund(ω₃|ω₁) and P_Exc(ω₃|ω₁) surfaces assists in connecting the spectroscopy to the potential shape. Qualitatively, as the potential asymmetry increases and the confinement decreases, the potential becomes more Eigen-like and the frequency ratio γ can take on values less than one. The extreme asymmetry narrows the potential around the zero point and first excited states, shifting ω₁₀ to greater than 1700 cm⁻¹. This means that the 2D IR spectrum at excitation frequencies below ω₁₀ = 1700 cm⁻¹ is primarily sampling the Zundel-like and intermediate sub-ensemble of configurations. Moreover, from Fig. 6, the Zundel-like vibrations appear more intense than the broad response from Eigen-like clusters. Additionally, the spectral response from intermediate geometric configurations resembles that of the Zundel-like clusters. Both of these factors suggest that IR spectroscopy more easily detects Zundel-like spectroscopic features over Eigen-like response, especially with the μ⁴ scaling in nonlinear experiments.

Comparing calculated 2D probability distributions with experimental 2D spectra indicates that the majority of the clearly resolved experimental features correlate with Zundel-like configurations and, to a lesser extent, intermediate configurations. However, in experiments, there are broad and featureless signatures in the spectral window 2000 cm⁻¹–2800 cm⁻¹, which are better explained in terms of Eigen-like configuration. These features have cross peaks to other resonances, which are consistent with Zundel-like and intermediate configurations. As a result, it is possible that fluctuations between a distribution of configurations including Eigen-like, intermediate, and Zundel-like configurations will best describe experimental 2D IR spectra. Given the lack of important dynamical effects in our calculations, and possible uncertainties in experimental spectra (such as water background subtraction), our conclusions must be viewed one step forward in the aqueous proton problem through a continuing advance of theoretical and experimental methods.

See the supplementary material for vibrational assignment tables of gas-phase and aqueous-phase clusters, additional probability density surfaces, histograms of aqueous-phase cluster geometries, additional one-dimensional and two-dimensional slices of the proton stretch PES, and additional comparisons between 2D IR spectra and conditional frequency probability surfaces.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors thank Rajib Biswas and Greg Voth for the MS-EVB trajectory. Q.Y. and J.M.B. thank the National Science Foundation, Grant No. CHE-1463552, for funding. W.B.C., J.H.H., and A.T. thank the Office of Basic Energy Sciences, U.S. Department of Energy, for support under Grant No. DE-SC0014305. B.D. thanks the Swiss National Science Foundation for support through the Postdoc.Mobility fellowship, Grant No. P400P2_180765.