Decoding the spectroscopic features and time scales of aqueous proton defects

The structure and transport of proton defects in aqueous solution underpin processes ranging from voltage-gated channels in biology to the development of improved proton exchange membrane materials.^1–4 While it has long been appreciated that proton defects form a wide range of structures whose interconversion enables rapid proton transport,^5–10 identification of the dominant motifs and their time scales remains a significant experimental^11–16 and theoretical^5,7,16–23 challenge.

However, of particular current interest is how these proton structures and their interconversion time scales manifest spectroscopically in the condensed phase.^{13,14,24–29} Linear vibrational spectroscopies provide one approach to probe these features. However, the electronic and structural flexibility of proton defects gives rise to a broad and largely featureless continuum between the far infra-red (IR) and the O–H stretch band regions of the absorption spectrum [see Fig. 1(a)]. Indeed, the flexibility of proton defects has been highlighted in gas-phase studies where even subtle chemical perturbations around the defect have been observed to induce significant spectral changes.^12,30,31 This raises the question as to whether the broad linear vibrational spectrum in the solution arises simply from rapid fluctuations of the proton defect structure or whether multiple defect structures exist that undergo slower interconversion with each contributing their own spectral signatures.

Separating these two pictures requires knowledge of the time scales for the interconversion of proton defect structures. Nonlinear spectroscopy provides a route to obtain this more detailed information by allowing the measurement of the correlation times of the frequencies associated with proton defects. Specifically, recent two-dimensional (2D) IR experiments on hydrochloric (HCl) acid solutions probed a stretch-bend cross-peak and observed a time scale of 480 fs, which was interpreted as a lower bound on the Zundel complex lifetime.¹⁴ While this result suggests that certain proton defect structures interconvert over time scales longer than previously observed,¹³ the interpretation of the cross-peak decay rested largely on assigning frequencies to the Zundel gas-phase structure within a cluster-centric paradigm.^{26–28,32–34} However, it remains unclear whether an interpretation of the liquid spectra based on the limiting Eigen and Zundel gas-phase cluster structures provides a complete picture of the structural motifs of condensed-phase proton defects and their interconversion time scales.

Here we use classical and path integral ab initio molecular dynamics (AIMD) simulations to decode the linear and nonlinear vibrational spectra of concentrated aqueous acid solutions. In particular, we show that one can define a simple coordinate that reveals how the vibrational spectrum and its evolution time scales arise from the asymmetry of the hydration structure of proton defects.

We performed classical and path integral AIMD simulations of concentrated HCl solutions and an aqueous excess proton in the NVT ensemble at T = 300 K under periodic boundary conditions using generalized gradient approximation (GGA) and hybrid density functional theory (DFT) to describe the interactions. Concentrated acid systems of concentrations 2M and 4M were simulated at their experimental densities.³⁵ For the excess proton system, a cubic box with sides of length 12.42 Å containing 64 water molecules (corresponding to a density of 999.3 kg m⁻³) and 1 additional proton was simulated, with the initial configurations obtained by protonating equilibrated bulk water. Initial structures for the 2M and 4M concentrated acid solutions were obtained from classical molecular dynamics simulations of aqueous chloride ions using the Amoeba2013 polarizable force field.³⁶ Two configurations spaced by 5 ns were extracted from the trajectories, the required number of water molecules was then protonated, and the positions of all protons were minimized. A total of 590 ps of classical simulations were performed for the 2M solution and 420 ps were performed for the 4M solution, composed of trajectories of ∼100 ps each. Simulations were performed using the i-PI program³⁷ using a multiple time scale (MTS) integrator of the reversible reference system propagator algorithm (r-RESPA) form.³⁸ The classical MTS simulations employed a 2.0 fs time step for the full forces and a 0.5 fs time step for the reference forces. Equilibration of 5 ps was performed for each trajectory using a local Langevin thermostat with a time constant of 25 fs, while production runs used a global stochastic velocity rescaling thermostat³⁹ with a time constant of 1 ps. Due to the global coupling of the latter thermostat, it causes a negligible perturbation to the system dynamics.⁴⁰ 

Path integral AIMD simulations were performed using ring polymer contraction (RPC)^41–43 with centroid contraction, i.e., contraction to P′ = 1 replicas using an MTS propagator with an outer time step of 2.0 fs and an inner time step of 0.25 fs. Simulations performed using contraction to P′ = 4 replicas gave results graphically indistinguishable from those using P′ = 1 for the probability distribution along the proton sharing coordinate (see Fig. S9 of the supplementary material). Thermostatted ring polymer molecular dynamics (TRPMD) path integral simulations^44–46 were performed by thermostatting the non-centroid normal modes using Langevin thermostats.⁴⁰ A total of over 130 ps of path integral simulations were performed using the revPBE0-D3 functional. The combination of the hybrid revPBE0-D3 functional with path integral simulations was employed due to the spurious effects that are encountered when quantizing GGA functionals such as revPBE-D3.⁴⁷ 

The full forces were evaluated using the CP2K program^48,49 at the DFT level of electronic structure theory using either the revPBE^50,51 GGA functional or the corresponding revPBE0^52,53 hybrid functional, with D3 dispersion corrections⁵⁴ added in both cases. Atomic cores were represented using the dual-space Goedecker-Teter-Hutter pseudopotentials.⁵⁵ Within the Gaussian plane wave (GPW) method,⁵⁶ Kohn-Sham orbitals were expanded in the TZV2P basis set, while an auxiliary plane-wave basis with a cutoff of 400 Ry was used to represent the density. Hybrid functional calculations employed a Coulomb operator truncated⁵⁷ at R_c = 6 Å and the auxiliary density matrix method⁵⁸ with the cpFIT3 fitting basis set. The self-consistent field cycle was converged to an electronic gradient tolerance of ϵ_SCF = 5 × 10⁻⁷ using the orbital transformation method⁵⁹ with the initial guess provided by the always-stable predictor-corrector extrapolation method^60,61 at each MD step. Full forces for the gas-phase Eigen and Zundel clusters were evaluated using the TeraChem program⁶² at the DFT level of electronic structure theory using the revPBE-D3^50,51 GGA functional for classical simulations and the revPBE0-D3^52,53 hybrid functional for path integral simulations.

The MTS and RPC reference forces were evaluated at the self consistent charge density functional tight binding (SCC-DFTB3)⁶³ level of theory using the DFTB+ program.⁶⁴ The 3ob parameter set was used for the H and O atoms⁶⁵ which was combined with a recently introduced parameterization for the hydrated halide ion.⁶⁶ Dispersion forces were included via a Lennard-Jones potential⁶⁷ with parameters taken from the Universal Force Field.⁶⁸ 

To identify the subensemble of protons connected to overcoordinated oxygen atoms, we employed the nearest-neighbor criterion, assigning each proton to the nearest oxygen atom. Overcoordinated oxygen atoms were then identified as those that were triply coordinated.

We first consider the experimental and simulated linear IR spectra of 4M HCl solution and neat water in Fig. 1(a). The spectra were calculated directly from the dipole correlation function obtained from ab initio molecular dynamics simulations. The quantum results were obtained from TRPMD simulations, which include quantum effects of the nuclei such as zero-point energy and tunneling, while the classical ones do not. While the classical simulations employed the GGA revPBE-D3 functional, due to the large zero-point energy contribution that is introduced when the nuclei are treated quantum mechanically (∼5 kcal/mol per hydrogen bond), the more accurate hybrid revPBE0-D3 functional was employed for the TRPMD simulation to avoid spurious effects that are encountered when quantizing GGA functionals.⁴⁷ Indeed, we have recently shown that, while classical simulations using the revPBE-D3 GGA functional appear to give the correct position of the O–H stretch band in water under ambient conditions⁴⁷ and at temperatures over the whole liquid range,⁶⁹ this arises from a fortuitous cancellation of error between the neglect of zero-point energy and an artificially weak chemical bond arising from its failure to correctly describe the electronic surface.⁴⁷ However, one should not rely on this cancellation in general, particularly in cases such as ours where the presence of proton defects and ions leads to a much wider range of hydrogen bond environments. Hence, to test this explicitly, we performed both classical simulations using the cheaper revPBE-D3 GGA functional and more advanced simulations combining path integrals with the hybrid revPBE0-D3 functional.

Both the simulated and experimental spectra exhibit continuous absorption between the water bend at ∼1600 cm⁻¹ and the beginning of the water O–H stretch at ∼3000 cm⁻¹. Taking the difference between the acid and water spectra allows one to identify distinct proton spectral features [Fig. 1(b)]. In particular, in addition to the increased intensity between the water bend and stretch, the bend broadens and a feature emerges around 1250 cm⁻¹. This feature, which has previously been suggested to be a Zundel-like shuttling motion,^11,12,70 corresponds to a proton motion orthogonal to the O–O axis in our simulations (see Fig. S1 of the supplementary material) and is thus associated with a bending motion of the proton defect, consistent with earlier simulation studies.^24,28 The ab initio simulations and experiments exhibit excellent agreement for both the total and difference spectra over the entire region for which the experimental data are available (>600 cm⁻¹), with the small discrepancy in the intensity depletion around 2250 cm⁻¹ arising from the underestimation of the bend-libration coupling feature in the simulated water spectrum. Indeed, for this particular observable, the previously discussed error cancellation between the GGA potential energy surface and nuclear quantum effects seems to hold remarkably well even for 4M acid solutions. The features of the IR spectrum, computed from the dipole autocorrelation function,⁷¹ arise from the underlying motions of the nuclei and are thus also encoded in the vibrational density of states (VDOS) of the system. This is shown in panels (c) and (d) of Fig. 1, where the absolute and difference VDOS spectra exhibit the same structural features as the IR, albeit with different intensities.

Given the correspondence between the VDOS and IR spectra, we therefore consider the time evolution of the VDOS of the protons in the acid solution to assess how proton spectral features arise and interconvert. As detailed in Sec. 1 of the supplementary material, we obtain the time-dependent spectrum as a sequence of instantaneous spectra, which are generated at each time by applying a smooth narrow window to the time series of the relevant quantity (velocity in the case of VDOS) and then performing a Fourier transform of that segment. Figure 2 shows the time evolution of the spectrum of one proton in a 4M HCl solution over a 100 ps simulation segment. During this segment, the proton is initially part of a water molecule (black), is part of a proton defect structure at ∼28 ps (orange), after which it returns to being part of a water molecule, and finally returns to being part of a proton defect (red). It then continues to fluctuate in and out of proton defect environments throughout the remainder of the simulation segment. The left-hand panel shows the spectrum associated with the proton at the three times indicated (black, orange, and red) compared to the total spectrum that proton experiences over the 100 ps simulation segment (gray). From this, one can see that the proton spectral features identified in Fig. 1 are observed when the proton forms part of a proton defect structure (orange and red). The time evolution of the VDOS thus provides a way to associate the vibrational spectrum of a specific proton over the entire frequency range and its evolution with the continually shifting local structural features in the condensed-phase environment of the acid solution.

The time evolution of the simulated spectrum can then be used to compute the correlation of vibrational intensities at different frequencies as a function of time in order to describe the transfer of vibrational energies between those two frequencies. Specifically, we calculate the correlation function corresponding to the pump (3150 cm⁻¹) and probe (1760 cm⁻¹) frequencies used in Ref. 14. Doing this for our aqueous excess proton defect simulation, one obtains a time constant of 1.4 ± 0.3 ps (see Fig. S3 of the supplementary material), which is consistent with the experimentally observed bound of >480 fs.¹⁴ The self-correlation time of the probe frequency (1760 cm⁻¹) obtained from our simulation is 1.6 ± 0.3 ps (see Fig. S4 of the supplementary material), suggesting that, consistent with the experimental interpretation, the primary mechanism for the slow relaxation time scale of that feature is via interconversion to a species that absorbs strongly at 3150 cm⁻¹.

What is the molecular origin of this time scale and what structural interconversion is it probing? To assess this, we consider the subensemble of protons directly connected to overcoordinated oxygen atoms (see Sec. II). There are three such protons by construction, which we refer to as the defect protons. While this allows us to define a subensemble of defect protons, it does not bias our subsequent analysis to any particular structure of the proton defect. Having done this, it is instructive to first consider the spectra of protons in defects where the hydrogen bond environment around the defect is explicitly asymmetric due to the formation of a hydrogen bond to a chloride ion (Fig. 3). Comparing the spectra of protons which are hydrogen bonded to the chloride (orange) and those that are not (blue), we observe that the former exhibit reduced, while the latter exhibit enhanced, proton spectral features, compared to the whole proton defect subensemble. This reflects the greater ability of the water oxygen atoms to share the defect protons, compared to the chloride ion, and thus underscores the impact of environment asymmetries on the spectral features of the proton defect. However, while the presence of a chloride ion in the proton defect’s 1st coordination shell explicitly breaks the symmetry of its chemical environment, the case of a defect coordinated solely by water is more subtle.

We now consider how broken symmetry around proton defects coordinated solely by water molecules manifests spectroscopically. In this case, all three acceptors are chemically identical in the 1st coordination shell of the defect (Fig. 4). One approach to differentiating the defect protons’ solvation environments is to consider the “hydrogen bond balance” of the water molecules they point at. In particular, the hydrogen bond balance is defined as the difference between the number of hydrogen bonds a water molecule receives and the number it donates (Fig. 4). This definition, which measures whether a water molecule is under or overcoordinated and hence indicates its propensity to accept or donate additional hydrogen bonds, has been successfully used to analyze the properties of numerous hydrogen bonded systems.^6,72–75 An illustration of this procedure is shown in Fig. 4 for three defect protons (blue, green, and orange), each of which points at a water with a different hydrogen bond balance arising from its (correspondingly colored) solvation environment. However, considering only the hydrogen bond balances of the acceptor water molecules at which defect protons point leads to poor separation of the spectrum. This is shown in Fig. S5 of the supplementary material, where the spectrum resolved only by the hydrogen bond balance gives near identical spectra in each case. The failure of this property alone to decompose the spectrum can be rationalized by considering cases where, although a defect proton might hydrogen bond to a water molecule with a high or low hydrogen bond balance, the other protons in the defect point at water molecules with the same hydrogen bond balances. One might expect that all these defect protons would exhibit similar sharing since none of them is directed at a water molecule that, based on its hydrogen bond balance, is a better acceptor than the others. The origin of the failure of the hydrogen bond balance to decompose the spectrum thus lies in its neglect of the relative chemical environments of the other defect protons.

Hence, we define a proton asymmetry coordinate, ϕ, that reports on the asymmetry of the hydrogen bond balances of water molecules in the 1st coordination shell of the defect. Specifically, for each defect proton, its asymmetry ϕ is calculated as the difference between the hydrogen bond balance of the water to which it is hydrogen bonded and the average hydrogen bond balance of the water molecules to which the other two defect protons are hydrogen bonded. For example, if all three defect protons point at water molecules with the same hydrogen bond balance, then this parameter yields a ϕ value of zero for each proton, indicating a lack of asymmetry for all of the protons. The calculation of the proton asymmetry coordinate for the protons in a defect is illustrated in Fig. 4. We emphasize that an asymmetry value, ϕ, is assigned to each defect proton and that the sum of the ϕ values for the three defect protons is zero by construction. It is also important that this coordinate does not assume a specific proton defect structure beyond the ability to identify overcoordinated oxygen atoms in the system. Although the illustration in Fig. 4 could be seen as resembling a solvated Eigen structure, since the ϕ parameter is a property of a particular proton and not the oxygen to which it is bound, it continues to be well defined even for protons that transfer. Therefore, for a proton that is highly shared between two oxygen atoms, one could equally well represent this figure as the extended Zundel-like depiction in Fig. S6 of the supplementary material. In what follows, we demonstrate that this purely structural property allows one to cleanly separate the proton spectral features and elucidate their interconversion time scales.

Figure 5 shows the defect proton vibrational spectrum resolved by the asymmetry coordinate, ϕ, for a simulation of an aqueous excess proton defect. The top panel shows the distribution of ϕ values obtained from the simulation, which is unimodal and centered around 0, with a range of values from −1.5 to +2.0. The spectrum for the lowest observed value of ϕ (yellow line in the left panel) exhibits large enhancements of the 1250 cm⁻¹ feature, bend broadening and density in the region between the water bend and stretch, which typify the proton vibrational spectral features observed in Fig. 1. This is because protons with very negative values of ϕ point at an acceptor water molecule that donates more hydrogen bonds than it accepts, whereas the other two protons of the defect are directed at water molecules with larger hydrogen bond balances. As such, a proton with a negative value of ϕ is hydrogen bonded to an excellent hydrogen bond acceptor water molecule, relative to the acceptors of the hydrogen bonds of the other defect protons. This results in preferential sharing of that proton with its first coordination shell water molecule, which manifests as enhanced proton spectral features. By contrast, the highest observed values of ϕ correspond to the opposite case, leading to a low preference for sharing. In this case, the spectrum (purple line in the left panel) exhibits few proton spectral features and indeed closely resembles that of the defect protons pointing at a chloride ion in Fig. 3. Hence, consideration of the relative asymmetry of all the hydrogen bonds formed by the defect, as encoded in ϕ, is essential to understanding the origin of the proton features of vibrational spectra.

In addition to allowing us to decode the linear spectrum, we can also use the structural parameter ϕ to rationalize the observed time scales observed in recent 2D-IR experiments.¹⁴ In particular, one can see from Fig. 5 that only protons with large values of ϕ (purple line) have strong spectral intensity at the pump frequency used in that experiment (3150 cm⁻¹), while the intensity is much greater at the probe frequency (1760 cm⁻¹) for low ϕ protons (yellow line). Hence, the 2D-IR experiment can be interpreted in terms of the ϕ coordinate as probing the time it takes for a proton to convert from a high- to low-ϕ structure, which is characterized by the ϕ correlation time. The ϕ correlation time extracted from our simulations (see Fig. S7 of the supplementary material) of an aqueous excess proton defect is 1.3 ± 0.3 ps, in excellent agreement with the frequency-frequency correlation time obtained from our simulations (1.4 ps) and consistent with the 2D-IR experiment (>480 fs). This time scale is similar to the orientational relaxation time of bulk liquid water observed experimentally (1.6-2.5 ps^76–79) and obtained from TRPMD simulations of revPBE0-D3 (1.6-2.6 ps depending on the molecular axis⁴⁷). The time scale extracted from the experiment can therefore be understood to arise from hydrogen bond rearrangements in the 1st and 2nd hydration shells of the proton defect that decorrelate ϕ by changing the overall asymmetry of its environment and which therefore gives rise to a time scale similar to the water reorientation time. The ϕ time correlation function, and hence its relaxation time scale, is essentially identical when it is calculated without allowing changes in the overcoordinated oxygen atom (see Fig. S7 of the supplementary material). This suggests that the primary mechanism for its decorrelation is water rearrangement around the proton defect rather than proton transfer, which leads to a change in the overcoordinated oxygen.^7,9,80 This analysis shows that the local structural properties of a hydrogen bond (e.g., the hydrogen bond balance) cannot provide a complete picture of the spectroscopy of aqueous protons, which is instead determined by the balance of hydrogen bonds around the defect.

We can now also assess how our approach relates to the commonly invoked paradigm of using the proton sharing coordinate δ, the difference of the distances of the proton from its two nearest oxygen atoms, to classify proton defect structures. Within this approach, structures with highly shared protons (δ close to zero) are commonly assigned as corresponding to “Zundel”-type structures, while other structures (δ further from zero) are assigned to “Eigen”-type or numerous other purported species.^{17,26,28,33,81} Figure 6 shows the distribution of δ resolved at each ϕ value for simulations with classical (top panel) and quantum (middle panel) nuclei. Although values of the proton sharing coordinate δ close to 0 are more prominent for negative ϕ values, there is a significant overlap between the distributions of δ, i.e., a given value of δ could arise when the proton is in many different solvation environments of the proton defect. In each of these environments, the proton rapidly explores a large range of δ values, which leads to the δ coordinate correlation function for the proton defect decaying to ∼0.2 in less than 100 fs (Fig. S8 of the supplementary material). This fast decay, corresponding to the exploration of the local environment of a particular proton defect structure, is consistent with that probed in recent ultrafast spectroscopy experiments.^13,82 After this massive and rapid decay, we observe a slower time scale with a time constant of 0.8 ± 0.1 ps, arising from the change in the environmental asymmetry as ϕ decorrelates (Fig. S7 of the supplementary material). Hence, the decay of the residual of the δ correlation function does subtly encode a slower time scale consistent with that probed in the 2D-IR experiments, although without providing any insight into the molecular rearrangements that lead to it.

It is also instructive to consider how the distribution along the proton sharing coordinate across the wide range of ϕ values in the liquid compares to those of the gas-phase Zundel and Eigen clusters at the same temperature, which are often invoked as the limiting proton defect structures that could exist in aqueous solution. The bottom panel of Fig. 6 shows the range of δ explored by the gas-phase Eigen and Zundel complexes at 300 K with classical and quantum nuclei. Comparing these to the ϕ-decomposed δ distributions in the top two panels, we can see that the highest proton asymmetry (ϕ) values visited in the acid solution have proton sharing (δ) distributions that resemble the Eigen cluster, while the lowest ϕ values show substantial population around δ = 0, as for the Zundel cluster. Indeed, it is remarkable that the hydrogen bond asymmetry, which arises from the differing numbers of hydrogen bonds received and donated from water molecules in the 2nd coordination shell of the proton defect, is enough to give rise to such a broad range of sharing patterns and spectral features of the defect.

Finally, although our results indicate that the particular pump-probe relaxation time probed in recent 2D-IR experiments¹⁴ is primarily associated with changes in the solvation environment around the proton defect, which gives rise to changes in proton sharing patterns, it is instructive to consider how a structural motif that has been associated with proton transfer events will manifest in the vibrational spectrum. Fig. S10 of the supplementary material shows the VDOS resolved by whether or not the proton defect is accepting a hydrogen bond. The formation of this hydrogen bond to the defect has been shown to be an important step in the proton transfer mechanism in previous simulations.⁹ While Fig. S10 of the supplementary material shows that the presence of the donated hydrogen bond to the defect does not change the vibrational spectrum of the proton defect’s protons, the H atom that is donated to the defect has a blue-shifted O–H stretch compared to both the defect protons and even those in liquid water. This hints that one may be able to probe the existence and time scales of the accepted hydrogen bond to proton defects by using pump and/or probe pulses in this high frequency region.

Employing our recent methodological advances,⁴³ which have enabled us to access nanosecond time scales, we have shown that ab initio molecular dynamics simulations of concentrated acid solutions exhibit excellent agreement with both the experimental IR spectra and 2D-IR relaxation time scales. By decomposing the time-dependent vibrational spectrum of the proton defects across the entire range of frequencies in their full condensed-phase chemical environment, we showed that both the proton spectral features and experimental spectroscopic time scales can be decoded using a structural parameter, ϕ, which encapsulates the overall solvation asymmetry of a proton’s chemical environment. Rather than attempting to rationalize the spectroscopy by considering the properties of individual hydrogen bonds, this parameter provides a collective description of the asymmetry of the proton’s solvation, which allows it to distinguish a broad range of proton sharing environments. Perhaps most importantly, it is the extrema of the asymmetry coordinate that give rise to the spectroscopic signatures that coincide with the pump-probe frequency combination studied in recent 2D-IR experiments. The >480 fs experimental time scale can therefore be interpreted as emerging from the relaxation of the proton asymmetry coordinate, which occurs by interconversion through a wide range of solvation structures that arise from water hydrogen bond rearrangements in the 2nd coordination shell. This provides a physically transparent framework for decoding the relationship between proton defect structures and their spectroscopic features and time scales in complex chemical environments.

See supplementary material for details of the calculation of the time-dependent VDOS, the identification of the 1250 cm⁻¹ spectral feature, validation of the width of the window used to obtain the VDOS, VDOS intensity cross- and self-correlation functions for the frequencies of interest, further illustration of the ϕ coordinate, time correlation functions of the δ and ϕ coordinates, and ring polymer contraction convergence validation.

We greatly thank William Carpenter and Andrei Tokmakoff for providing the experimental linear acid spectra and for helpful comments on the manuscript. This material is based upon work supported by the National Science Foundation under Grant No. CHE-1652960. T.E.M. also acknowledges support from a Cottrell Scholarship from the Research Corporation for Science Advancement and the Camille Dreyfus Teacher-Scholar Awards Program. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We would also like to thank Stanford University and the Stanford Research Computing Center for providing computational resources and support that have contributed to these research results. This work used the XStream computational resource, supported by the National Science Foundation Major Research Instrumentation program (No. ACI-1429830).