Nature of Excess Hydrated Proton at the Water–Air Interface

Understanding the interfacial molecular structure of acidic aqueous solutions is important in the context of, e.g., atmospheric chemistry, biophysics, and electrochemistry. The hydration of the interfacial proton is necessarily different from that in the bulk, given the lower effective density of water at the interface, but has not yet been elucidated. Here, using surface-specific vibrational spectroscopy, we probe the response of interfacial protons at the water–air interface and reveal the interfacial proton continuum. Combined with spectral calculations based on ab initio molecular dynamics simulations, the proton at the water–air interface is shown to be well-hydrated, despite the limited availability of hydration water, with both Eigen and Zundel structures coexisting at the interface. Notwithstanding the interfacial hydrated proton exhibiting bulk-like structures, a substantial interfacial stabilization by −1.3 ± 0.2 kcal/mol is observed experimentally, in good agreement with our free energy calculations. The surface propensity of the proton can be attributed to the interaction between the hydrated proton and its counterion.

(a) shows the off-resonant SFG response from the H2O-air interface for pure H2O and H2O containing 1M HCl, NaCl, NaI, and NaOH in the subphase as a function of infrared frequency in the OH stretch vibrational region. Indeed, neat H2O does not show any resonant contribution. The frequency-independent trend of the SFG response can be found even in the 1M NaCl, NaI, and NaOH solutions, while the 1M HCl solution shows a frequency- Frequency (cm -1 ) S3 dependent SFG feature that is an extended part of the resonance around 3000 cm -1 . This means that the non-resonant response does not change at 1M ionic strength and the changes in the hydrogen-bonded OH stretch region in the presence of HCl must be an effect of presence of protons at the surface. Figure S1 (b) shows the SFG responses from the D2O-air interface containing 1M DCl, NaCl and NaOD in the subphase as a function of infrared frequency in the OH-rather than ODstretch vibration region. Similar to H2O, D2O also does not show any resonant contribution. Moreover, the non-resonant SFG response does not change in the presence of 1M DCl, NaCl and NaOD. This means that the non-resonant response in presence of protons at the interface is the same as that of pure water, suggesting that the rise of the SFG response in the presence of 1M HCl may be an effect of the appearance of new vibrational modes and/or a change of dispersion. Figure S2: The SFG intensity in SSP polarization from the H2O-air interface in the presence of 1M NaCl in the subphase.

Fitting procedure
We quantify the amplitudes of the different resonances in the SFG spectra using an established fitting procedure. According to this fitting procedure, the ISFG is proportional to the square of the second-order susceptibility χ (2) of the sample and intensity of the incoming infrared (IR) and visible (VIS) beams, The ( ) is a sum of a non-resonant terms, ( ) (comprising of a non-resonant amplitude and a non-resonant phase ) and resonant ( ( ) ) contribution(s). Each ( ) is expressed in terms of a Lorentzian line shape with amplitude , central frequency , and half-width at half maximum .
The phase of the interfacial proton-continuum The broad featureless SFG feature throughout the range ~2000-3000 cm -1 increases in the presence of 1M HCl in the subphase (Fig 1, main text), similar to the enhancement of IR and Raman response for HCl in the bulk water. Following the studies in the bulk, this enhancement is identified as an enhancement due to the presence of different vibrational modes of the Eigen and Zundel moieties of the proton. The enhancement in the SFG intensity response can be an associated modulation in the imaginary and/or real part of the susceptibility (Im ( ) and Re ( ) ). Indeed, while the intensity of the imaginary part of the susceptibility seems to be finite only for frequencies above ~2900 cm -1 , in the intensity spectrum "proton continuum" seems to span a much wider frequency range. To address this apparent contradiction, we compared intensity spectra with the real and imaginary responses. Figure S3: (a-c) Imaginary and real part of χ (2) spectra as a function of frequency in SSP polarization from pure H2O-air interface and in the presence of 1M HCl in the subphase. (d-f) Experimentally measured SFG intensity overlaid with the intensity spectra constructed from the experimental imaginary and real χ (2) spectra Figure S3 (a-c) represent the imaginary (solid lines) and real (dotted lines) part of the ( ) responses as a function of frequency in SSP polarization from the water-air interface of pure H2O and H2O containing HCl with 1M concentration in the subphase. Within the experimental noise level, the Im ( ) response of H2O seems to be nearly zero in the range ~2350-2800 cm -1 and clearly negative in the range ~2800 -3500 cm -1 . The real part of the response in the presence of 1M HCl is definitely more negative than that of pure H2O in the range ~2350-2900 cm -1 , implying that the enhanced SFG intensity in this region is dominated by the Re ( ) contribution. This implies that the extended intensity observed in the intensity spectrum originates from the real part of the susceptibility.

Figures 3 (d-f
) represent the experimental SFG response overlaid with the constructed intensity from the experimental ( ) spectra. The constructed intensity spectra match fairly well with the experimental spectra bolstering the fidelity of the data.

Angular distribution of the orientation of the free OD:
As shown in figure 3 in the main text, the intensity of the free OH/ OD decreases in the presence of acid in the subphase. This can be a result of a reduction in the number of free OH/OD oscillators, but also the result of a change in the average angular distribution -or a combination of the two. The , respectively) in the following ways: where ( = , , )is the Fresnel factor and is the incidence or reflection angle of the light of frequency with respect to the surface normal along the z-axis, to the surface in the x-y S6 plane. The amplitudes of the free OH/OD peak and in the ( ) and ( ) spectra, respectively, are related to 〈 〉/〈 〉 via in the slow motion limit. 1,2 For the slow motion limit, the decay of the orientational memory of the free OH/OD is much slower than its vibrational relaxation. r = 0.15 is given by the ratio of the transition polarizabilities ̇ ∥ ⁄ for the free OH stretch mode, 3 where ∥ and ⊥ denote directions parallel and perpendicular to the OH bond, respectively.
To check whether the angular distribution changes with varying the proton concentration, we have performed SFG measurement with SSP a n d PPP polarization combinations at the D2O-air interface for pure D2O and D2O containing 2M DCl in the subphase. The intensities are then fitted with an established Lorentzian model to extract the areas of the free OD vibration.  Figure S4(a-b) shows SFG intensities at the neat D2O-air and 2M DCl D2O solution-air interfaces in SSP and PPP polarization combinations, respectively. The spectra shown with red circles correspond to the experimental intensity spectra. The black lines show the fit for the experimental spectra. Figure 4(c-f) show the imaginary and real parts of ( ) (in black) calculated from the fits shown in panel (a) and (b).
Each spectrum has a sharp peak at ~2720 cm -1 and a broad response that extends into the low-frequency region. The sharp feature is the vibration of the non-hydrogen bonded OD oscillators that point away from the surface into the air (free OD group). The broad feature is the vibrational response of hydrogen-bonded OD groups that broadly vary in hydrogen-bond strength. Compared to that of pure D2O, the intensity of the free OD peak decreases and the hydrogen-bonded OD stretch peak increases upon the addition of 2M DCl into the subphase.
Since, the experimental intensity spectra comprise of interfering resonant vibrations of different signs, we fitted the experimental homodyne spectra using an established procedure under the constraint that the shape of the imaginary and real parts of ( ) obtained from the fits matches the experimental spectra. The areas of the free OD stretch peaks were used to calculate the ratios, / . The calculated ratios are 0.42 and 0.41 for D2O and 2M DCl solution, respectively.  Figure S5 (a) shows the calculated area of free OD stretch vibration in SSP and PPP polarizations as a function of the average angle that free OD makes with the surface normal. It shows that the area in SSP polarization depends on the angles very differently from that in PPP. Figure S5 (b) shows the calculated / from the / as a function of average angle of the free OD with respect to the surface normal in slow dynamics limit. The ratios for D2O and 2M DCl (0.42 and 0.41 respectively) correspond to the angles 61.3 0 and 57.5 0 , i.e. a difference of ~4 0 which is negligible given the noise of the measurement and uncertainties of the fits. Therefore, our results show that the average angular distribution of free OD is unchanged upon the presence of acid, and the variation of the SFG signal can be thus attributed solely to the variation of the number of free OD stretch chromophores. Figure S6: (a-c) SFG spectra and fits with the Lorentzian lineshape model in the OD stretching range for D2O solutions with varying NaCl and HCl concentrations for three different sets of experiments of which the average is shown in Fig. 3 of the main text. The total ion concentration is the same for all measurements: 1M.

Computational protocol for ab initio molecular dynamics simulation
Ab initio molecular dynamics (AIMD) simulations were carried out by using the CP2K software package 4 based on Born-Oppenheimer propagation via the Quickstep electronic structure module. The Quickstep approach combines a Gaussian basis set for the atomic orbitals with an auxiliary plane wave basis for the electron density. We chose a triple-ζ quality TZV2P basis set and use a charge density cutoff of 320 Ry together with the NN50 smoothing for the charge density and spline smoothing for the derivative of the charge density. We employed the revised PBE functional (revPBE) 5 combined with Grimme's empirical dispersion D3(0) correction. 6 The core electrons were represented by the norm-conserving Goedecker-Teter-Hutter pseudopotentials. 7 The choice of these calculation methods arises from the recent evaluation of the description of the interfacial water where we concluded that the revBPE-D3(0) nicely reproduces the properties of bulk and surface water among the generalized gradient approximation level of theory. 8

S9
We simulated the slab HCl(aq) system containing 10 H + and Clion pairs together with 160 water molecules. This provides a concentration of 3.4 M, higher than the experimental condition of 1M. We used this high concentration to highlight the impact of the H3O + and Clions on the interfacial water structure. We used the simulation cell of 16.3 Å × 16.3 Å × 44.1 Å and the periodic boundary conditions. The z-axis forms the surface normal and the xy-plane is parallel to the surface. The AIMD simulations are performed in the canonical (NVT) ensemble using the CSVR thermostat, 9 where the target temperature was set to 300 K. We integrated the equation of motion by using a time step of 0.5 fs. 10 initial configurations are prepared, and subsequently they were equilibrated for 10 ps. After 10 ps equilibrations, we ran 60 ps AIMD simulations for sampling the trajectory. Thus, we obtained a total of 10 × 60 ps trajectories, which were used for analyzing the structure of the interface and compute the SFG signal.

Definition of interface
Density profiles of each atomic/molecular species as well as the whole system along the z-axis are shown in Figure S7. It is clearly shown that H + (aq) is surface-active whereas Cl -(aq) is inactive, as consistent with force field molecular dynamics simulations. To define the position of interface, the density profile of the whole system is fitted by the hyperbolic tangent, where is the z-coordinate of the target atom, is the z-coordinate of the Gibbs dividing surface, δ is defined as the interfacial thickness parameter, and denotes the bulk density. We also used the instantaneous liquid interface description, 10 where the liquid interface is determined based on the individual snapshots. The comparison of the density profiles of H3O + ions between the average description and instantaneous liquid interface description is shown in Figure S8.

Computational protocols for SFG spectra
Although several computational protocols to calculate SFG responses of water (or aqueous solutions)-air interfaces have been proposed, almost all of these methods are not applicable or practically infeasible for the protonated water-air interface. Subsequently, we calculated the SFG responses of H + (aq) at water-air interfaces based on the normal mode calculations of protonated water clusters, instead of calculating the SFG response directly from the time correlation function approach using AIMD trajectories. 11 We selected the fully capped H + (aq) cluster, (H5O2 + ) 4H2O, at water-air interface from these AIMD trajectories and calculated their SFG responses by using quantum chemistry calculations. The sampling and computational procedures are as follows: (1) We sampled the molecular configurations from the AIMD trajectories every 1 ps. (2) We assigned each hydrogen atom to the oxygen atoms which are closest to the hydrogen atom. In this manner, we defined 150 H2O and 10 H3O + molecules. (3) We calculated the δ-coordinate, which is defined by the difference in distances between the intra-and inter-molecular O-H distances δ = |OH1 -OH2|, for all O-H bonds of H3O + . A small (large) δ-coordinate means that the proton is shared by the two oxygen atoms with similar (biased) weight. That means that the proton with the smallest δ-coordinate is the most active for each H3O + . When the z-coordinate of the most active proton satisfies | − | < 3.11 Å, a H3O + ion is considered to stay at water-air interfaces. (4) For the most active O-H bond, we further picked up its H-bond acceptor water molecule from the AIMD configuration. The selected water molecule and the H3O + molecule constitute the H5O2 + species. (5) We further capped four hydrogen atoms of non-active O-H bonds according to the AIMD snapshot. In this manner, we sampled multiple conformations of the Eigen type (H3O + , δ~0.5) and the Zundel type (H5O2 + , δ~0.0) almost in similar local environments.

S11
Once we sampled the (H5O2 + ) 4H2O configurations at the water-air interface, we relaxed coordinates of only hydrogens/proton of H5O2 + and performed normal mode calculations to evaluate the SFG responses of H + (aq) at the water-air interface. ssp-SFG response in normal mode approximation is expressed by 12 Here, and or denote the z-component of dipole moment and the xx(yy)-element of the polarizability of the whole cluster, respectively, and Qi is the i-th normal mode. Totally more than 2500 clusters were sampled and SFG responses were calculated. Although the coordinate of (H5O2 + ) ⋅4H2O is taken from AIMD trajectories, only H5O2 + keeps the original Hbonds and other four water molecules loose H-bonds except from H5O2 + . Thus, we consider only normal modes in which the contribution from H5O2 + exceeds half. The resultant SFG responses of O-H stretch in H5O2 + were binned with a frequency resolution of 100 cm -1 to plot spectra. All normal mode calculations were performed by using revPBE-D3/aug-cc-pVDZ level theory with the ORCA quantum chemistry package. 13 In the AIMD simulation at the revPBE-D3(0)/TZV2P level of theory, the ratio of the Eigen and Zundel structures is 3:7. However, the Eigen structure is known to be underestimated in the simulation with classical nuclei. Thus, we corrected the number of the Eigen and Zundel moieties by multiplying the factor of ∆ / , where ∆ is the potential barrier artificially obtained from the simulation with classical nuclei. This ∆ was estimated to be 0.6 kcal/mol. 14 With this correction, we estimated that the ratio of the Eigen and Zundel structure is 1:1. This correction factor is also used for spectral calculation.

Sensitivity check of the Eigen/Zundel definition to SFG spectra
Since the transition from the Eigen structure to the Zundel structure is contentious, there is no unique criteria of the δ-coordinate to disentangle these two species in the simulated H 3O + species. In the main text, we employed the criteria of δ > 0.3 for Eigen and δ < 0.1 for Zundel according to Ref. 15. In Figure S9, we also plot the SFG responses for Eigen and Zundel moieties where the criterion of Eigen and Zundel is set at δ > 0.2 and δ < 0.2, respectively. Shapes of these spectra are basically identical to those in the main text and, thus, the spectra of the Eigen and Zundel structures is insensitive to the threshold.   16 From the AIMD trajectories, we obtained the angle formed by the free O-H group and the surface normal of 66 degrees for the HCl solution-air interface. Since the free O-D group angle at D2O-air interfaces was also 66 degrees, 17 the data indicates that the invariance of the angle due to the presence of HCl. This is consistent with the orientation analysis based on the SFG free O-D data in the SSP and PPP polarization combinations.

Simulated Free O-H group
Furthermore, we computed the fraction of the free O-H/O-D group at these interfaces. We observed the free O-H/D-D group with the surface density of 1.77 nm -2 and 2.51 nm -2 for the DCl solution-air and D2O-air interfaces, respectively. 17 These are now compared with the experimental data. This is plotted in Figure S10. The experimental and simulation data show good agreement, indicating that the simulation captures the trend observed in the experiment accurately.

Hydration structure of H3O + ion
To examine the hydration structure of H3O + ion at the water-air interface and in the bulk, we computed the RDF of the H3O + ion and the hydrogen atom of water molecule. The data is shown in Figure S11. This indicates that the neighboring hydration structure is very similar between the water-air interface and the bulk. Only oxygen atoms which are in the range of |rz -rG| < 3.11 Å are considered for surface, whereas oxygen atoms in the rage of |rz| < 3.11 Å are used for bulk.

Proton dynamics at interfaces
We computed the accommodation time of the proton at the interfacial region (|rz -rG| < 3.11 Å) via; where P(t) = 1 when the oxygen atom of the H3O + ion at time t is in the interfacial region, 0 otherwise. The data is displayed in Figure S12. This shows a lifetime of ~13 ps. This is consistent with the previous study. 18 Figure S13 displays the time evolution of the individual H3O+ ions. This indicates that hydrated protons in the bulk region diffuse much faster than at interface.

S15
Additional simulation data Figure S14. The simulated depth profile of water density within the instantaneous liquid interface description and the averaged interface description.